曲線描画

(3) B-スプライン曲線

ベジェ曲線は、バーンスタイン基底関数をブレンディング関数として使ったパラメトリック曲線であることを前回紹介しました。ブレンディング関数は他にも様々な形式を考えることができますが、今回はその中の一つとして「B-スプライン基底関数(B-Spline Basis Function)」をブレンディング関数とした「B-スプライン曲線(B-Spline Curve)」を紹介します。

1) B-スプライン基底関数 (B-Spline Basis Function)

「B-スプライン基底関数 (B-Spline Basis Function)」は、「De Boor-Coxの漸化式」と呼ばれる以下の漸化式で表される関数です。

但し、二番目の式において、N_j,p(t) = N_j^p と表しています。t_j は「ノット (Knot)」と呼ばれる数列で、ここでは要素数を M + 1 とします。すなわち、0 を開始番号とすれば、{ t_j } = ( t₀, t₁, ... t_M ) となります。但し、i < j ならば t_i ≤ t_j である必要があります。また、p はこの基底関数の次数を表します。この式のなかで、ノット列の添字の最大は j + p + 1 なので、M ≥ j + p + 1 が成り立ち、基底関数は N_0,p(t) から N_M-p-1,p(t) までの M - p 個になります。

p = 0 のときは、上に示した漸化式は使えません。この時は、次のように定義されます。

N_j,0(t) による線形結合はいわゆる「階段関数 (Step Function)」と呼ばれる不連続な関数です。

図 1-1. Nj0 のグラフ (階段関数)

M は p より大きければどんな値でも取りえます。また、p は M に依存することなく決めることができます。p = 1 の場合、漸化式から

となるので、N_j,1(t) は t = t_j のとき 0 となり、線形に増加しながら t = t_j+1 で 1 に、その後今度は線形に減少して t = t_j+2 で再び 0 になります。よって N_j,1(t) のグラフは、底辺が t_j から t_j+2 で、高さが 1 の三角形になります。

図 1-2. Nj1 のグラフ (三角形関数)

p = 2 になると、式がかなり複雑になってきます。そこで、t_j+m = t_m と表して j は書かないようにして、さらに

と表すと、

と三つの区分に分かれた二次式で表現されます。t = t_m のとき u_m,n = 0、t = t_n のとき u_m,n = 1 になることに注意すると、t = t₀, t₃ のとき N_j,2(t) = 0 であり、t_m - t_n = D_m,n で表すと、t = t₁ のときは

になります。N_j,2( t₁ ) の計算には、t₀ ≤ t ≤ t₁, t₁ ≤ t ≤ t₂ のどちらの定義域の式を使っても結果が同じであることに注意して下さい。同様に、t = t₂ のときは

になります。定義域 t_m ≤ t ≤ t_m+1 での N_j,2( t ) を N_j,2[m]( t ) で表すと、N_j,2[0]( t ), N_j,2[2]( t ) はそれぞれ

なので、t = t₀, t₃ で極値を持つことがわかります。従って (二次項の係数は正なので)、t₀ ≤ t ≤ t₁ では短調増加、t₂ ≤ t ≤ t₃ では単調減少であり、もし N_j,2( t ) が連続なら t₁ ≤ t ≤ t₂ では極値を持たなければなりません。実際、この定義域での極値を求めると、

より、

なので、

と計算できて、左辺がゼロの時は

になります。この値と t₁ の差は

であり、t₂ との差は

なので、確かに極値は必ず t₁ ≤ t ≤ t₂ の範囲内に存在します。また、

と表せることを利用して、この導関数に t = t₁, t₂ をそれぞれ代入すると、

となり、N_j,2[0](t), N_j,2[2](t) の導関数が

であることから、各式に t = t₁, t₂ をそれぞれ代入すると、N'_j,2[0]( t₁ ) = N'_j,2[1]( t₁ ), N'_j,2[1]( t₂ ) = N'_j,2[2]( t₂ ) であることが容易にわかります。従って、N_j,2(t) の導関数は t = t₁, t₂ にて連続であることになります。しかし、二階導関数は各区分において定数となり、

なので、連続にはなりません。

図 1-3. Nj2 のグラフ

p = 3 の場合、

となります。ここでも、t = t_m のとき u_m,n = 0、t = t_n のとき u_m,n = 1 になることを利用すれば、

であり、しかも各区分の境界で連続であることが容易にわかります。また、

となるので、それぞれの式は三重根 t₀, t₄ を持ち、この点は変曲点になります。N_j,3[0]( t ) は三次項の係数が正なのに対し、N_j,3[3]( t ) は負であることから、N_j,3[0]( t ) は t ≥ t₀ で短調増加、N_j,3[3]( t ) は t ≤ t₄ で短調減少です。さらに、u_m,n・u_r,s・u_v,w の導関数が

であることを利用すると、N'_j,3[1]( t ) は

となるので、再び u_m,n( t_m ) = 0, u_m,n( t_n ) = 1 を利用して N'_j,3[1]( t₁ ), N'_j,3[1]( t₂ ) を計算すると

となります。同様に、N'_j,3[2]( t ) が

であることを利用して N'_j,3[2]( t₂ ), N'_j,3[2]( t₃ ) を計算すると、

という結果が得られ、N'_j,3[1]( t₂ ) = N'_j,3[2]( t₂ ) であることが証明できます。また、

なので、p = 2 の場合と同様に、導関数も t = t₁, t₂, t₃ において連続であることになります。

さらに二階導関数を求めてみましょう。N⁽²⁾_j,3[0]( t ), N⁽²⁾_j,3[3]( t ) は、

なので、

です。N⁽²⁾_j,3[1]( t ), N⁽²⁾_j,3[2]( t ) は、u_m,n・u_r,s・u_v,w の二階導関数が

であることを利用して

となるので、

と求めることができて、t = t₁, t₂, t₃ では二階導関数も連続であることを示すことができます。三階導関数は定数で、

なので、これは連続にはなりません。

図 1-4. Nj3 のグラフ

かなり苦労して p = 3 までの関数の様子を確認しました。この結果を見ると、以下のことが推察できます。

N_j,p( t ) の漸化式は次のようなものでした。

N_j,p( t ) は、N_j,p-1( t ) と N_j+1,p-1( t ) の和によって計算されます。N_j,0( t ) は t_j ≤ t ≤ t_j+1 の範囲のみにゼロ以外の正値 1 を持った階段関数であり、N_j,p( t ) は最終的に N_j,0( t ) で表すことができるので、ノットによって多項式が区分けされるのは明らかです。また、p = 1 において、t ≤ t_j, t_j+2 ≤ t の領域では N_j,1( t ) = 0 でした。つまり、t_j から t_j+2 以外の領域において N_j,1( t ) が値を持たないことを意味します。N_j,2( t ) が、N_j,1( t ) と N_j+1,1( t ) の和で表されるということは、前者が t_j から t_j+2、後者が t_j+1 から t_j+3 までの領域のみ値を持っていることから、t_j から t_j+3 までの領域のみ値を持つということになります。これを繰り返し考えれば、N_j,p( t ) が t_j から t_j+p+1 までの領域のみ値を持つということが示されます。

t = t_j のとき、漸化式の右辺第一項は明らかにゼロです。また第二項も、N_j+1,p-1( t_j ) = 0 が成り立てばゼロになります。同様に、t = t_j+p+1 のとき、漸化式の右辺第二項はゼロになり、第一項は、N_j,p-1( t_j+p+1 ) = 0 が成り立てばゼロになります。N_j+1,p-1( t ) は t_j+1 から t_{(j+1)+(p-1)+1} までの領域以外はゼロであり、N_j,p-1( t ) は t_j から t_j+(p-1)+1 までの領域以外はゼロなので、N_j+1,p-1( t_j ) = N_j,p-1( t_j+p+1 ) = 0 であり、N_j,p( t_j ) = N_j,p( t_j+p+1 ) = 0 であることが証明できます。

p = 1 のとき、任意の t に対して N_j,1( t ) ≥ 0 になります。p > 2 の任意の p に対して N_j,p-1( t ) ≥ 0 が成り立つと仮定したとき、漸化式において右辺第一項にある ( t - t_j ) / ( t_j+p - t_j ) は t ≥ t_j においてゼロ以上であり、第二項にある ( t_j+p+1 - t ) / ( t_j+p+1 - t_j+1 ) は t ≤ t_j+p+1 においてゼロ以上なので、t_j ≤ t ≤ t_j+p+1 の範囲で N_j,p( t ) ≥ 0 が成り立ちます。t_j ≤ t ≤ t_j+p+1 以外の範囲では N_j,p-1( t ) も N_j+1,p-1( t ) もともにゼロですから、全ての領域において N_j,p( t ) ≥ 0 が成り立ちます。従って、帰納法により任意の p に対して N_j,p( t ) はゼロ以上の値をとることになります。

N_j,p-1( t ) と N_j+1,p-1( t ) が連続であると仮定した時、[ ( t - t_j ) / ( t_j+p - t_j ) ]・N_j,p-1( t ) と [ ( t_j+p+1 - t ) / ( t_j+p+1 - t_j+1 ) ]・N_j+1,p-1( t ) はやはり連続であり、その和である N_j,p( t ) も連続です。N_j,1( t ) が連続なので、これも帰納法から p ≥ 1 の任意の p に対して N_j,p( t ) は連続であることが示されます。N_j,p( t ) が連続であるのは、正値をとる範囲 t_j ≤ t ≤ t_j+p+1 のみではないことに注意して下さい。N_j,p( t_j ) = N_j,p( t_j+p+1 ) = 0 が必ず成り立つことから、全領域に対して N_j,p( t ) は連続です。

さらに N_j,p( t ) の p - 1 階までの導関数が連続であることを証明するためには、一階導関数に対する以下の漸化式が成り立つことを示す必要があります。

N'_j^p+1 は、積の導関数の公式を利用して

と展開することができますが、先ほどの漸化式が成り立つと仮定したとき、

とさらに展開することができます。項の数が増えたので再び t_j+m を t_m で表し、(1), (2) の和に対する N_j^p-1 と N_j+1^p-1 の係数をそれぞれ K₁, K₂ としてその値を求めると

(3), (4) の和に対する N_j+1^p-1 と N_j+2^p-1 の係数をそれぞれ K₃, K₄ とすると

になり、K₂ + K₃ は

となることから、N'_j^p+1 は最終的に

となり、p で漸化式が成り立てば、p + 1 でも同じ式が成り立つことになります。p = 1 のとき、導関数 N'_j,1( t ) は、

なので漸化式が成り立ち、帰納法により、ゼロより大きな任意の p に対しても成り立つことになります。

導関数に対する漸化式から、N'_j,p( t ) は N_j,p-1( t ) と N_j+1,p-1( t ) の線形結合で表されることがわかります。p > 1 ならば、N_j,p-1( t ) と N_j+1,p-1( t ) は連続です。また、N_j,p-1( t ) は [ t_j, t_j+p ] のみ値を持ち、N_j,p-1( t_j ) = N_j,p-1( t_j+p ) = 0 が成り立ちます。同様に、N_j+1,p-1( t ) は [ t_j+1, t_j+p+1 ] のみ値を持ち、N_j+1,p-1( t_j+1 ) = N_j+1,p-1( t_j+p+1 ) = 0 なので、それらの線形結合である N'_j,p( t ) は連続で、[ t_j, t_j+p+1 ] のみ値を持ち、N'_j,p( t_j ) = N'_j,p( t_j+p+1 ) = 0 が成り立ちます。

図 1-4. N'j3 と Nj2, Nj+12 のグラフ

二階導関数は N'_j,p-1( t ) と N'_j+1,p-1( t ) の線形結合で表されることになるので、上記と同様に考えれば、二階導関数 N⁽²⁾_j,p( t ) は連続であり、[ t_j, t_j+p+1 ] のみ値を持ち、N⁽²⁾_j,p( t_j ) = N⁽²⁾_j,p( t_j+p+1 ) = 0 であることがわかります。これを繰り返すと、0 ≤ k ≤ p - 1 の範囲において k 階導関数は連続であり、[ t_j, t_j+p+1 ] のみ値を持ち、N^(k)_j,p( t_j ) = N^(k)_j,p( t_j+p+1 ) = 0 であることが示されます。しかし、p 階導関数では成り立たなくなります。例えば p = 1 のとき、上式をみれば明らかなように N'_j,1( t ) は不連続です。N⁽²⁾_j,2( t ) は N'_j,1( t ) と N'_j+1,1( t ) の線形結合ですから、これも連続になるとは限りません。

N_j,p( t ) は「基底関数」という名が付けられていることから推測できるように、基底として利用することができます。N_j,p( t ) による線形結合を

としたとき、P_M,p( t ) = 0 が任意の t で成り立つためには全ての係数 a_j がゼロになる必要があることが示されれば N_j,p( t ) は基底関数になります。p = 0 のとき、

より、任意の t に対して P_M,p( t ) がゼロになるためには明らかに全ての係数 a_j がゼロでなければなりません。N_j,p-1( t ) が基底関数であると仮定したとき、

より

と表されるので、

が任意の t に対して成り立たなければなりません。最初の式によって、a₀ = 0 である必要があることは明らかです。j = 1 の場合、第二式に a₀ = 0 を代入すると a₁( t - t₁ ) = 0 が任意の t に対して成り立つ必要があることになって、a₁ = 0 でなければなりません。これを繰り返すと、任意の j に対して a_j = 0 である必要があることになり、N_j,p( t ) も基底関数になります。よって、帰納法によって、任意の p に対して N_j,p( t ) は基底関数です。

ところで、a_j が全て 1 だったとき、1 ≤ j ≤ M - p - 1 における N_j^p-1 の係数は

となるので、

となります。N_0,p-1( t ) が値を持つのは t₀ ≤ t ≤ t_p の範囲であり、N_{M-(p-1)-1,p-1}( t ) が値を持つのは t_M-p ≤ t ≤ t_M の範囲なので、t_p ≤ t ≤ t_M-p の範囲においては N₀^p-1 と N_{M-(p-1)-1,p-1}( t ) の項は無視することができて

となります。実際に、p = 1 の場合を求めてみると、

となりますが、N_j⁰ は t_j ≤ t ≤ t_j+1 のみ 1 で残りの範囲はゼロなので、

という結果が得られます。これは、t₁ ≤ t ≤ t_M-1 において Σ_j{0→M-2}( N_j¹ ) = 1 であることを意味し、p = 2 に対して同じ処理を行うと t₂ ≤ t ≤ t_M-2 の範囲で P_M,2( t ) = Σ_j{1→M-3}( N_j¹ ) となるので、N₀¹ が t₀ ≤ t ≤ t₂、N_M-2¹ が t_M-2 ≤ t ≤ t_M の範囲でしか値を持たないことから

になります。これを繰り返し適用することによって、任意の p に対して

という結果を得ることができます。

B-スプライン基底関数を求めるためのサンプル・プログラムを以下に示します。

/*
  CalcBSplineBasisFunc : B-スプライン基底関数の計算

  const vector<double>& knot : ノット列
  unsigned int j : ノット列の開始番号
  unsigned int p : 次数
  double t : 計算対象の独立変数

  ノット列は昇順である必要があるが、そのチェックは行わない

  戻り値 : 計算結果
*/
double CalcBSplineBasisFunc( const vector<double>& knot, unsigned int j, unsigned int p, double t )
{
  if ( knot.size() == 0 ) return( NAN );

  // ノット列のデータ長が充分でない場合は nan を返す
  unsigned int m = knot.size() - 1;
  if ( m < j + p + 1 ) return( NAN );

  // 正値をとる範囲外ならゼロを返す
  if ( ( t < knot[j] ) || ( t > knot[j + p + 1] ) ) return( 0 );
  // p = 0 かつ knot[j] <= t <= knot[j + p + 1] なら 1 を返す
  if ( p == 0 ) return( 1 );
  // p = 1 の場合、三角の頂点の値は特別扱い
  if ( p == 1 && t == knot[j + 1]  ) return( 1 );

  // 漸化式の計算
  double d1 = ( knot[j + p] == knot[j] ) ? 0 :
    ( t - knot[j] ) * CalcBSplineBasisFunc( knot, j, p - 1, t ) / ( knot[j + p] - knot[j] );
  double d2 = ( knot[j + p + 1] == knot[j + 1] ) ? 0 :
    ( knot[j + p + 1] - t ) * CalcBSplineBasisFunc( knot, j + 1, p - 1, t ) / ( knot[j + p + 1] - knot[j + 1] );

  return( d1 + d2 );
}

内容は非常にシンプルで、「deBoor-Coxの漸化式」を使って次数の低い式から順に再帰的に計算をしているだけです。p = 0 なら階段関数になるので、値のある領域内なら 1 を返します。また、p の値に関係なく、t が値を持つ領域を外れている場合は即座にゼロを返しています。

なお、N_j,0( t_j+1 ) = N_j+1,0( t_j+1 ) = 1 にすると、漸化式による N_j,1( t_j+1 ) の計算結果は

となってしまい、本来の値 1 とは異なる結果となります。それを防ぐため、N_j,1( t_j+1 ) の計算は特別扱いして常に 1 を返すようにしています。

漸化式を計算するときの注意点として、t_j = t_j+p となるようなノットが含まれた場合、漸化式の項の中でゼロ除算が発生します。このときは、t_j ≤ t_j+1 ≤ ... ≤ t_j+p-1 ≤ t_j+p が成り立たなければならないことから t_j = t_j+1 = ... = t_j+p-1 = t_j+p であることを意味します。t_j = t_j+p のとき、次数が p の B-スプライン基底関数でゼロ除算が発生するのは

のときで、N_j-1,p( t ) のときは第二項、N_j,p( t ) のときは第一項がゼロ除算を含みます。それらは N_j^p-1 を含みますが、その定義域は t_j ≤ t ≤ t_j+p であり、存在しないものとみなしても問題はないことになります。よって、ゼロ除算が発生する項は無視する(ゼロとする)対処が必要です。

2) B-スプライン曲線 (B-Spline Curve)

B-スプライン基底関数をブレンディング関数とした曲線を「B-スプライン曲線 (B-Spline Curve)」といいます。M 個の点列 { P_j } ( j = 0, 1, ... M - 1 ) に対し、B-スプライン曲線 S(t) は以下の式で表されます。

N_M-1,p( t ) は、t_M-1 ≤ t ≤ t_(M-1)+p+1 のみの領域で正値をもつので、ブレンディング関数の数は M なのに対し、必要なノット列は t₀ から t_M+p までの M + p + 1 個になることに注意して下さい。また、t₀ ≤ t ≤ t₁ においては S( t ) の和の成分は N_0,p( t ) のみであり、N_0,p( t₀ ) = 0 であることから、S( t₀ ) = 0 となります。また、同様の考え方によって S( t_M+p ) = 0 です。したがって、t₀ ≤ t ≤ t_M+p の範囲で曲線を描画した場合は、必ず原点でつながった閉曲線の形になります。

後ほど説明をしますが、ベジェ曲線のブレンディング関数であるバーンスタイン基底関数は B-スプライン基底関数の特殊形です。従って、ベジェ曲線の持つ特性の多くは B-スプライン曲線にもあります。

B-スプライン基底関数が基底となることは前節ですでに証明しています。また、t_p ≤ t ≤ t_M の範囲においては Σ_j{0→M-p-1}( N_j,p( t ) ) = 1 であり、N_j,p( t ) がゼロ以上の値であることもわかっているので、t_p ≤ t ≤ t_M の範囲で描画すればアフィン不変性と凸閉包性が成り立ちます。

対称性に対しては、

が成り立てば

となるので、点列を P_M-1 から P₀ へ、t を t_M+p から t₀ へと逆側にたどっても曲線は同じになります。しかし、ベジェ曲線については対称性が常に成り立つのに対し、B-スプライン曲線の場合は常に成り立つとは限りません。例えば p = 0 のとき、N_M-1-j,0( t₀ + t_M - t ) は t₀ + t_M - t が t_M-1-j から t_M-j の間にあるときだけ 1 になるので、

より

で 1 になります。N_j,0( t ) が t_j から t_j+1 の間で 1 となることから、

でなければならず、ノット列を両端から中央方向へ順番に抽出して求めた和が全て等しくなる必要があります。例えば、ノット列を中央で二つに分けた時、前半と後半でそれぞれ同値となったり、ノット列が自然数列などの等差数列になればこの条件を満たします。

t = t_p のとき S( t_p ) = P₀ となるためには、N_0,p( t_p ) = 1、N_j,p( t_p ) = 0 ( j > 0 ) である必要があります。しかし、これは明らかに、常には成り立ちません。同様に、S( t_M ) = P_M-1 も常には成り立たないので、B-スプライン曲線は一般には端点を通過しません。しかし、p = 1 のときは常に N_j,1( t_j+1 ) = 1 であり、N_0,1( t₁ ) = N_M-1,1( t_M ) = 1 なので、必ず端点を通過するのみでなく、全ての点を通過することも保証されます。すなわち、p = 1 の場合の B-スプライン曲線は各点を結んだ折れ線になります。

以上の結果だけでは、ベジェ曲線に比べて特性としては劣っているようにみえますが、ベジェ曲線にはない優れた特性として「局所制御」があります。ベジェ曲線の場合、点列の一つを移動すると曲線全体にその影響が及びます。ベジェ曲線では基底関数が全定義域にわたって値を持ち、係数の変化が全定義域に及ぶためにそうなるのですが、B-スプライン基底関数の定義域は次数 p で決まるので、ある点 P_k が変化した時に影響を受ける基底関数は t_k を定義域に含むもの、すなわち N_k-p-1^p から N_k^p までとなります。しかし、N_k-p-1^p と N_k^p は t_k が定義域の端にあたり、そこでの値はゼロなので、曲線の軌跡には影響を与えません。結局、影響のある基底関数は N_k-p^p から N_k-1^p までの p 個となります。また、S( t ) が変化する定義域は t_k-p ≤ t ≤ t_k+p となります。
ベジェ曲線では新たな点を追加した上でいくつかの区分に分けてスプライン曲線として描画をすることで局所制御できるように対応していましたが、B-スプライン曲線の場合は p を制御することによって局所性を持たせながら曲線を描画することができます。ベジェ曲線の場合、端点で接線となるため、スプライン曲線として描画するときに分割する二点(一つのスプライン曲線の終点と、次のスプライン曲線の始点のこと)を結んだ線分上に新たな点を追加してそれを各スプライン曲線の終点・始点とすれば、二つのスプライン曲線は滑らかに接続されます。これは、それぞれの端点における接線の方向が等しくなることから成り立つのですが、さらに高階の導関数については必ずしも成り立ちません。実際、ベジェ曲線のブレンディング関数であるバーンスタイン基底関数の導関数は

なので、二階導関数は

となって、t = 0 のとき

t = 1 のとき

なので、ベジェ曲線の t = 0, 1 での二階導関数の値は

となります。但し、P_m,n は、点 P_m から点 P_n までを結んだベクトルを表します。点 P_n と 点 P_n+1 を結んだ線分上に新たな点 P_m を用意し、P_m を終点及び始点としてスプライン曲線の形でベジェ曲線を描画する場合、二つのベジェ曲線は P_m で滑らかにつながるのでした。しかし、P_m における両曲線の二階導関数の値は N( N - 1 )( P_n,n-1 + P_n,m ) と N( N - 1 )( P_n+1,m + P_n+1,n+2 ) なので、P_n,n-1 + P_n,m = P_n+1,m + P_n+1,n+2 が成り立たない限り一致しません。

図 2-1. P_m における二階導関数の向きと大きさ

曲線において、k 階導関数の値が連続である場合は「C^k-連続 (C^k-Continuity)」といいます。曲線が全てつながっていれば C⁰-連続であり、接線の向きと大きさが連続であれば C¹-連続となります。上述のようにベジェ曲線をスプライン曲線として描画した場合、曲線は滑らかになりますが、スプライン曲線どうしがつながった点 P_m での接線ベクトルは N( P_m - P_n ) と N( P_n+1 - P_m ) なので、P_m - P_n = P_n+1 - P_m が成り立たなければ、方向が等しいだけで連続にはならず、C¹-連続にはなりません。このように、k 階導関数ベクトルの方向だけが等しく、大きさは異なるような場合は「G^k-連続 (G^k-Continuity)」といいます。従って、ベジェ曲線によるスプライン曲線は G¹-連続ということになります。ベジェ曲線では二階導関数ベクトルが向きも大きさも一致するとは限らないことから、C²-連続でも G²-連続でもありません。C^k-連続や G^k-連続は曲線の滑らかさを表す指標となり、大きな k に対して成り立つほど滑らかさの度合いは高くなります (ベジェ曲線のブレンディング関数であるバーンスタイン基底関数は高階導関数が存在するので、点列全てを使ってベジェ曲線を描画した場合は充分滑らかな曲線が得られます。上記の説明は、あくまでも点列をいくつかのグループに分割してスプライン曲線の形でベジェ曲線を描画した場合であることに注意して下さい)。

B-スプライン基底関数 N_j,p( t ) は、任意の t に対し、0 ≤ k ≤ p - 1 の整数 k について k 階導関数まで連続でした。従って、N_j,p( t ) の線形結合である S( t ) も k 階導関数まで連続であることになり、B-スプライン曲線は C^p-1-連続なので、ベジェ曲線を使ったスプライン曲線よりも滑らかな曲線を描画することができます。

B-スプライン曲線は、ノットの配列の仕方によって三つのタイプに分類されます。

1. 一様 (Uniform)

ノット列が等間隔に配置された状態の曲線

2. 開一様 (Open Uniform)

ノット列の両端が、p + 1 分だけ同じ値である場合の曲線

3. 非一様 (Non-uniform)

一様でも開一様でもないノット列の曲線

ノット列の間隔が等しければ、漸化式の一部は定数化することができます。例えば、t_j = j ( j = 0, 1, ... M ) である場合、

となります。p = 0 のときは幅が一定の矩形になり、形状は全て等しくなります。p = 1 ならば

であり、全ての j に対して底辺が 2 で高さ 1 の三角形となります。一様のときは、任意の p に対し、B-スプライン基底関数は定義域が異なるだけで全て同じ形状になります。従って、N_0,p( t ) さえ求めてしまえば、形状をシフトするだけで任意の j に対する基底関数が得られます。

図 2-2. 一様 B-スプライン基底関数 ( p = 2 )

開一様は、ノット列の両端の値が p + 1 だけ等しい場合を指します。開一様なノット列の最初の値を a とした時、t₀ = t₁ = ... = t_p = a なので、N_0,p( t ) は

となり、N_p⁰ の定義域が [ t_p, t_p+1 ] = [ a, t_p+1 ] なので N_0,p( t_p ) = N_0,p( a ) = 1 になります。同様に、ノット列の最後の値を b とした時の N_M-1,p( t ) は、t_M = t_M+1 = ... = t_M+p = b より

となり、N_M-1⁰ の定義域が [ t_M-1, t_M ] = [ t_M-1, b ] なので N_M-1,p( t_M ) = N_M-1,p( b ) = 1 になります。この結果は、開一様である場合は端点を通過するということを表します。ここで、ノット列の数を 2p + 2 とし、前半の p + 1 個を 0、残り後半を 1 としたとき、p = 0 ならばノット列は { 0, 1 } であり

となります。p = 1 ならば、ノット列は { 0, 0, 1, 1 } で

であり、任意の p に対しては t₀ = t₁ = ... = t_p = 0、t_p+1 = t_p+2 = ... = t_2p+1 = 1 なので

という結果が得られます。バーンスタイン基底関数の漸化式は

なので、N_j,p( t ) と b_j+1,p( t ) は 1 ≤ j ≤ p - 1 の間で右辺が一致します。左辺は添字が一つ分シフトしているので、下図 2-3 を見ると N_0,p( t ) と等しくなるバーンスタイン基底関数は b_p,p( t ) となるようにみえますが、実際には

より、N_0,p( t ) と b_0,p( t ) が一致します。同様に、N₁⁰ から N_p+1⁰ までの式で再帰的に得られる N_1,p( t ) と ( 実際には存在しませんが ) b₁⁰ から b_p+1⁰ までの式で再帰的に得られる b_p+1,p( t ) が一致し、b_p+1,p( t ) は b_1,p( t ) を意味します。さらに、N_p⁰ から N_2p⁰ までの式で再帰的に得られる N_p,p( t ) と b_p⁰ から b_2p⁰ までの式で再帰的に得られる b_2p,p( t ) が一致し、b_2p,p( t ) は b_p,p( t ) を意味します。よって、この場合の N_j,p( t ) はバーンスタイン基底関数そのものであり、描画される曲線もベジェ曲線そのものとなります。

図 2-3. B-スプライン基底関数とバーンスタイン基底関数の再帰的表現

下のグラフは、ノット列を { 0, 0, 0, 1, 1, 1 } としたときの二次の B-スプライン基底関数を表しています。

図 2-4. 開一様 B-スプライン基底関数 ( p = 2 )

これは、二次のバーンスタイン基底関数

と同一であることがこのグラフからもわかります。

非一様なノット列を利用すると、表現できる曲線の形状は増えます。例えば、二次の B-スプライン関数に対し、ノット列の途中に、値が等しいノットを次数と等しい二個だけ挿入すると、その開始番号を j としたとき、t_j = t_j+1 なので、t_m - t_n = D_m,n と表すと

となり、N'_j-1,2[1]( t ) は定義できません。このとき、導関数は t_j で不連続で、正値 1 / ( t_j - t_j-1 ) から負値 1 / ( t_j - t_j+2 ) に変わります。一般に、次数 p の B-スプライン曲線に対してノット列の途中に p 個の同数を挿入すると「尖点(Cusp)」として表現することができるので、急激に折れ曲がり尖った部分を作りたい時に有効です。

図 2-5. 非一様 B-スプライン基底関数 ( p = 2 )

上に示したグラフの左側は、ノット列を { 0, 0.2, 0.5, 1.5, 2.2, 3.0 } としたときの二次の B-スプライン関数です、右側は、中央の値を少し変えて { 0, 1.0, 1.5, 1.5, 2.2, 3.0 } としたもので、j = 1 のときのグラフは最大値をとるところで線が折れ曲がった状態になっています。この性質を利用すると、曲線に尖点を持たせることができるわけです。

B-スプライン曲線の軌跡を計算するためのサンプル・プログラムを以下に示します。

/**********************************************************
  BSplineCurveBase : B-スプライン曲線描画用基底クラス
**********************************************************/
class BSplineCurveBase : public ParametricEquation
{
  unsigned int p_;   // 次数
  vector<double> x_; // 点列の x 座標
  vector<double> y_; // 点列の y 座標

 protected:

  vector<double> knot_; // ノット列(protected)

 private:

  // calc : 曲線の X/Y 座標値を計算する(純粋仮想関数)
  virtual double calc( const vector<double>& p, double t ) = 0;

  // set : メンバ変数の初期化
  void set( vector<double>::const_iterator knot0,
            vector<double>::const_iterator knot1,
            vector< Coord<double> >::const_iterator p0,
            vector< Coord<double> >::const_iterator p1,
            unsigned int p );

 public:

  /*
    BSplineCurveBase コンストラクタ
  */

  /*
    配列の範囲を指定

    vector<double>::const_iterator knot0, knot1 : ノット列の範囲
    vector< Coord<double> >::const_iterator p0, p1 : 点列の範囲
    (knot1,p1 は末尾の次を表すため含まれない)
    unsigned int p : 次数
  */
  BSplineCurveBase( vector<double>::const_iterator knot0,
                    vector<double>::const_iterator knot1,
                    vector< Coord<double> >::const_iterator p0,
                    vector< Coord<double> >::const_iterator p1,
                    unsigned int p )
  { set( knot0, knot1, p0, p1, p ); }

  /*
    配列全体を指定

    const vector<double>& knot : ノット列
    const vector< Coord<double> >& pnt : 点列
    unsigned int p : 次数
  */
  BSplineCurveBase( const vector<double>& knot,
                    const vector< Coord<double> >& pnt,
                    unsigned int p )
  { set( knot.begin(), knot.end(), pnt.begin(), pnt.end(), p ); }

  // 点列のサイズを返す
  unsigned int size() const
  { return( x_.size() ); }

  // 次数を返す
  unsigned int order() const
  { return( p_ ); }

  // x, y の値を求める
  virtual double x( double t )
  { return( calc( x_, t ) ); }
  virtual double y( double t )
  { return( calc( y_, t ) ); }

  // 仮想デストラクタ
  virtual ~BSplineCurveBase() {}
};

/*
  BSplineCurveBase::set : メンバ変数の初期化

  vector<double>::const_iterator knot0, knot1 : ノット列の範囲
  vector< Coord<double> >::const_iterator p0, p1 : 点列の範囲
  (knot1,p1 は末尾の次を表すため含まれない)
  unsigned int p : 次数
*/
void BSplineCurveBase::set( vector<double>::const_iterator knot0,
                            vector<double>::const_iterator knot1,
                            vector< Coord<double> >::const_iterator p0,
                            vector< Coord<double> >::const_iterator p1,
                            unsigned int p )
{
  p_ = p;

  // 点列の大きさ
  unsigned int pntSize = std::distance( p0, p1 );
  // ノット列の大きさ
  unsigned int knotSize = std::distance( knot0, knot1 );

  // ノット列の初期化
  // 用意したノット列の要素が pntSize + p_ + 1 より少なければ、末尾側は最後の要素で補う
  // 逆に多ければ末尾側の要素は無視する
  if ( knotSize < pntSize + p_ + 1 ) {
    knot_.assign( knot0, knot1 );
    for ( unsigned int i = knotSize ; i < pntSize + p_ + 1 ; ++i )
      knot_.push_back( knot_.back() );
  } else {
    knot_.assign( knot0, knot0 + pntSize + p_ + 1 );
  }

  // 点列を X/Y に分解して初期化
  while ( p0 != p1 ) {
    x_.push_back( p0->x );
    y_.push_back( p0->y );
    ++p0;
  }
}

/**********************************************************
  BSplineCurve : B-スプライン曲線描画用クラス
                 B-スプライン基底関数をそのまま利用
**********************************************************/
class BSplineCurve : public BSplineCurveBase
{
  // calc : 曲線の X/Y 座標値を計算する
  virtual double calc( const vector<double>& p, double t );

 public:

  /*
    BSplineCurve コンストラクタ

    BSplineCurveBase をそのまま利用
  */

  BSplineCurve( vector<double>::const_iterator knot0,
                vector<double>::const_iterator knot1,
                vector< Coord<double> >::const_iterator p0,
                vector< Coord<double> >::const_iterator p1,
                unsigned int p )
    : BSplineCurveBase( knot0, knot1, p0, p1, p ) {}

  BSplineCurve( const vector<double>& knot,
                const vector< Coord<double> >& pnt,
                unsigned int p )
    : BSplineCurveBase( knot, pnt, p ) {}
};

/*
  BSplineCurve::calc : 曲線の X/Y 座標値を計算する

  const vector<double>& pnt : 点列の X/Y 座標
  double t : 計算対象の独立変数

  返り値 : 計算結果
*/
double BSplineCurve::calc( const vector<double>& pnt, double t )
{
  double ans = 0; // 計算結果
  unsigned int p = order(); // 次数

  for ( unsigned int j = 0 ; j < size() ; ++j ) {
    double d = pnt[j] * CalcBSplineBasisFunc( knot_, j, p, t );
    if ( ! isnan( d ) )
      ans += d;
  }

  return( ans );
}

ベジェ曲線描画のサンプル・プログラムの場合と同様に、B-スプライン曲線描画でも基底クラス BSplineCurveBase を用意し、純粋仮想関数 calc によって曲線の軌跡の計算方法が切り替えられるようにしています。メンバ変数は、次数 p_ と点列の X, Y 座標 x_, y_、ノット列 knot_ の 4 つあり、最後の knot_ は限定公開メンバです。
BSplineCurveBase のメンバ関数 set はメンバ変数を初期化するためのもので、コンストラクタから呼び出して利用します。点列のサイズ pntSize に対し、ノット列は pntSize + p_ + 1 以上のサイズである必要があるので、不足している場合はノット列の末尾の要素を使って補います。逆に多すぎる場合は末尾側の要素は無視して、ノット列のサイズがちょうど pntSize + p_ + 1 に等しくなるように調節しています。

BSplineCurve は、B-スプライン基底関数の計算処理関数 CalcBSplineBasisFunc を利用して曲線の軌跡を求めるためのクラスで、BSplineCurveBase からの派生クラスになっています。仮想関数の calc は、CalcBSplineBasisFunc を利用してB-スプライン基底関数の値を求め、点列の x, y 座標を係数として多項式を求めています。

次数 p、点列の数 M の B-スプライン曲線の軌跡は通常 t_p ≤ t ≤ t_M の範囲のみを描画します。それ以外の範囲を含めた場合、アフィン不変性と凸閉包性が成り立たなくなり、特に全てのノット列を含めた定義域で描画すると、曲線は必ず原点を含む閉曲線となってしまいます。このクラスは ( CalcBSplineBasisFunc 自身も ) 描画する範囲は限定していないため、t としてどのような値を渡しても何らかの値を返します。特に M < p となった場合、アフィン不変性と凸閉包性が成り立つ曲線は描画できなくなります。

渦巻状に並べた九つの点列に対し、一様ノット列を使って曲線を描画した結果を以下に示します。ノット列は 0 から始まり、次数 p に対して p + 9 までの、p + 10 個の要素からなる自然数列になります。また、描画する範囲は、最後の曲線を除いて、アフィン不変性と凸閉包性が成り立つ t_p ≤ t ≤ t₉ までとしています。
p = 1 の場合は各点を結んだ線分になり、端点だけではなく全ての点列が直線上に含まれますが、p > 1 の場合、曲線は端点を通りません。最後の描画結果は、p = 3 の場合に対してノット列の全領域 t₀ ≤ t ≤ t₁₂ で軌跡を求めたもので、左上原点で曲線の両端がつながった形になります。

図 2-6. 一様 B-スプライン曲線

p = 1	p = 2	p = 3	p = 3 (t0 ≤ t ≤ t12)

同じ点列を開一様ノット列で描画した結果は以下のようになります。p + 10 個のノット列に対し、前半と後半の p + 1 個の要素はそれぞれ 0 と 2 で、中間の値を 1 としています ( 実際のノット列の値は各図形の下側に示してあります )。

図 2-7. 開一様 B-スプライン曲線

p = 3	p = 4	p = 8	ベジェ曲線
0 0 0 0 1 1 1 1 1 2 2 2 2	0 0 0 0 0 1 1 1 1 2 2 2 2 2	0 0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2

p = 3 の場合、曲線は途中で途切れた形になります。j = 4 のとき、N₄³ はゼロ以外の値を持つ範囲が全くなく無視することができるので、B-スプライン基底関数は j = 0, 1, 2, 3 と j = 5, 6, 7, 8 が左右対称となるような形状になり、しかも t = 1 のところで完全に分断された形になります。

図 2-8. ノット列 { 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2 } の B-スプライン基底関数 ( p = 3 )

上図のグラフから明らかに、t = 1 のところでは N_3,3( 1 ) = N_5,3( 1 ) = 1 であり、S( 1 ) は P₃・N₃³ か P₅・N₅³、つまり P₃ と P₅ そのものと等しくなり、ここで曲線は切断された形になります。p = 4 のときは切断はされていませんが、同じ値が次数と等しい 4 つ分だけ並んでいるために尖点が発生しています。p = 8 ではノット列に 9 個の 0 と 2 を並べており、このときはベジェ曲線と同じ形状になっています。

下図は、ノット列をランダムに作成して非一様二次 B-スプライン曲線として描画した例です。

図 2-9. 非一様二次 B-スプライン曲線

(1)	(2)	p = 8
0.98 1.46 2.10 3.08 3.40 4.21 4.89 5.73 6.26 6.26 7.16 8.04	0.80 1.10 1.77 2.05 2.05 2.52 3.02 3.40 3.40 4.18 4.62 5.55	0.61 0.77 1.23 1.88 2.76 2.79 3.31 3.46 4.27 5.04 5.32 6.19

B-スプライン基底関数の漸化式

を使い、バーンスタイン基底関数に対する「ド・カステリョのアルゴリズム」導出のときと同様に、点列 P_j を係数とする B-スプライン多項式 S(t) を変形すると以下のようになります。

定義域を [ t_p, t_p+1 ] に限定すると、S(t) の値に影響する B-スプライン基底関数は N₀^p から N_p^p までのみになります。そこで、この定義域に限定すると、p = 1 の場合は

となるので、S(t) は線分 P₀ - P₁ を t - t₁ : t₂ - t の比率で分けたときの内点(重心)に相当します。さらに p = 2 の場合、

より、S(t) は 線分 P₀ - P₁ の t - t₁ : t₃ - t の内点と P₁ - P₂ の t - t₂ : t₄ - t の内点を結んだ線分をさらに t - t₂ : t₃ - t の比率で分けた内点になります。p = 3 の場合、

となるので、S(t) は 線分 P₀ - P₁ の t - t₁ : t₄ - t の内点と P₁ - P₂ の t - t₂ : t₅ - t の内点を結んだ線分を t - t₂ : t₄ - t の比率で分けた内点、P₁ - P₂ の t - t₂ : t₅ - t の内点と P₂ - P₃ の t - t₃ : t₆ - t の内点を結んだ線分を t - t₃ : t₅ - t の比率で分けた内点を結び、その線分を t - t₃ : t₄ - t の比率で分けた内点に相当します。この結果から、任意の p に対する S(t) を [ t_p, t_p+1 ] の範囲について求めるためには、まず P_j-1 と P_j ( j = 1, ... p ) を結んだ線分を t - t_j : t_j+p - t の比率で分けた p 個の内点 P_j-1⁽¹⁾ を求め、次に P_j-1⁽¹⁾ と P_j⁽¹⁾ ( j = 1, ... p - 1 ) を結んだ線分を t - t_j+1 : t_j+p - t の比率で分けた p - 1 個の内点 P_j-1⁽²⁾ を求める、という具合に、1 個の内点 P₀^(p) になるまで繰り返し求めればよいことがわかります。[ t_p+1, t_p+2 ] の定義域に対しては点列を一つずらして P₁ から P_p+1 までを使えばよく、最後の [ t_M-1, t_M ] に対しては P_M-p-1 から P_M-1 までを使って求めます。このようにして曲線の軌跡を求めるアルゴリズムを「De Boor-Cox のアルゴリズム (De Boor-Cox Algorithm)」といいます。ノット列を p + 1 個のゼロと p + 1 個の 1 にした場合は、1 ≤ j ≤ p に対して t - t_j : t_j+p - t = t : 1 - t が成り立ち、これは「ド・カステリョのアルゴリズム (De Casteljau's Algorithm)」と同じ操作になります。すなわち、「ド・カステリョのアルゴリズム」は「De Boor-Cox のアルゴリズム」の特殊形であると考えることができます。

図 2-10. De Boor-Cox のアルゴリズム

「De Boor-Cox のアルゴリズム」を使って B-スプライン曲線の軌跡を計算するサンプル・プログラムを以下に示します。

/******************************************************************
  DeBoorCox : De Boor-Cox のアルゴリズムを使った B-スプライン曲線
******************************************************************/
class DeBoorCox : public BSplineCurveBase
{
  vector<double> buff_; // 内点の計算結果を保持するバッファ

  // calc : 曲線の X/Y 座標値を計算する
  virtual double calc( const vector<double>& p, double t );

 public:

  /*
    DeBoorCox コンストラクタ

    BSplineCurveBase をそのまま利用
  */

  DeBoorCox( vector<double>::const_iterator knot0,
             vector<double>::const_iterator knot1,
             vector< Coord<double> >::const_iterator p0,
             vector< Coord<double> >::const_iterator p1,
             unsigned int p )
    : BSplineCurveBase( knot0, knot1, p0, p1, p ),
      buff_( p ) {}

  DeBoorCox( const vector<double>& knot,
             const vector< Coord<double> >& pnt,
             unsigned int p )
    : BSplineCurveBase( knot, pnt, p ),
      buff_( p ) {}
};

/*
  DeBoorCox::calc : 曲線の X/Y 座標値を計算する

  const vector<double>& pnt : 点列の X/Y 座標
  double t : 計算対象の独立変数

  返り値 : 計算結果
*/
double DeBoorCox::calc( const vector<double>& pnt, double t )
{
  if ( size() <= 1 ) return( NAN );
  unsigned int p = order(); // 次数

  // t を含むノット列の領域 [ knot_[s-1], knot_[s] ] の添字 s を求める
  unsigned int s = 0;
  for ( ; s < knot_.size() ; ++s )
    if ( t < knot_[s] ) break;
  if ( s <= p || s > size() )
    return( NAN );
  s -= p;

  // p = 0 のときは点列の座標値そのものを返す
  if ( p == 0 ) return( pnt[s - 1] );

  // 点列から最初の内点を求めて buff_ に格納する
  vector<double>::iterator it = buff_.begin();
  for ( unsigned int j = s ; j < p + s ; ++j )
    *it++ = ( ( knot_[j + p] - t ) * pnt[j - 1] + ( t - knot_[j] ) * pnt[j] )
      / ( knot_[j + p] - knot_[j] );

  // 再帰的に内点を計算
  for ( unsigned int e = p ; e > 1 ; --e ) {
    unsigned int inc = p - e + 1; // 前側のノットの添字に対する増分
    it = buff_.begin();
    for ( unsigned int j = s ; j < s + e - 1 ; ++j ) {
      *it = ( ( knot_[j + p] - t ) * *it + ( t - knot_[j + inc] ) * *( it + 1 ) )
        / ( knot_[j + p] - knot_[j + inc] );
      ++it;
    }
  }

  return( buff_[0] );
}

曲線の軌跡を求める前に、ノット列の中で t を含む領域の添字を求める処理を行います。得られた添字 s は、領域の終端を表すことに注意して下さい。また、このとき t が t_p から t_M ( M は点列のサイズ ) の範囲内にあるかをチェックし、そうでなければ nan を返しています。BSplineCurve クラスの場合とは異なり、「De Boor-Cox のアルゴリズム」では t_p ≤ t ≤ t_M 以外の領域の計算はできないため、計算処理前に必ずチェックをする必要があります。チェック後は p を減算し、点列の添字に合わせた値に切り替えます。

p = 0 のとき、計算式の分母はゼロになってしまうため、点列の座標をそのまま返します。そうでなければ、先ほど示したアルゴリズムに従って再帰的に内点の計算を繰り返し、内点が一つになったところでその値を返します。

一様・開一様・非一様のそれぞれに対し、三次の B-スプライン曲線を、BSplineCurve と DeBoorCox それぞれのクラスを使って描画した結果を以下に示します。ノット列は図の下側に示してあります。どちらの描画方法を使っても、全ての場合において曲線は一致します。

図 2-11. BSplineCurve と DeBoorCox それぞれの描画結果

	一様	開一様	非一様
BSplineCurve
DeBoorCox
ノット列	0 1 2 3 4 5 6 7 8 9 10 11 12	0 0 0 0 1 1 1 1 1 2 2 2 2	0 0.5 1.5 3.0 4.5 6 6 6 6.25 6.5 7.5 10 15

3) NURBS ( Non-Uniform Rational B-Spline )

B-スプライン曲線においては、ノット列をコントロールすることによって、曲線を任意の箇所で折れ曲がった状態にするなど、曲線を自由に描くことができることを前節で説明しました。ノット列は、一つ一つの B-スプライン基底関数の幅や形状を変化させることができます。

ここで、各基底関数に対する重み値 w_j を定義してスプライン曲線を表現してみます。

w_j = 0 ならば、w_j・P_j・N_j,p( t ) の項は存在しないことになり、P_j は曲線の形状に影響を与えないことになります。逆に、w_j の値が大きくなるほど P_j の影響力は大きくなります。

さて、このままだと Σ_j{0→M-1}( w_j・N_j,p( t ) ) ≠ 1 なので、アフィン不変性が成り立たなくなります。そこで、全体を Σ_j{0→M-1}( w_j・N_j,p( t ) ) で割ることで正規化します。よって、スプライン曲線は

となります。これを「非一様有理 B-スプライン ( Non-Uniform Rational B-Spline ; NURBS )」といいます。

B-スプライン曲線の描画ができるのなら、NURBS の実装は非常に簡単です。早速サンプル・プログラムを示したいと思います。

/**************************************
  NurbsCurve : NURBS 曲線描画用クラス
**************************************/
class NurbsCurve : public BSplineCurveBase
{
  std::vector< double > w_; // 重み値

  // calc : 曲線の X/Y 座標値を計算する
  virtual double calc( const std::vector< double >& p, double t );

  // setWeight : 重み値の取り込み
  void setWeight( std::vector< double >::const_iterator w0,
                  std::vector< double >::const_iterator w1 );

 public:

  /*
    NurbsCurve コンストラクタ
  */

  /// 配列の範囲を指定して構築
  ///
  /// knot0,knot1 : ノット列の範囲 (knot1 は末尾の次を表すため含まれない)
  /// p0,p1 : 点列の範囲 (p1 は末尾の次を表すため含まれない)
  /// w0,w1 : 重み値の範囲 (w1 は末尾の次を表すため含まれない)
  /// p : 次数
  NurbsCurve( std::vector< double >::const_iterator knot0,
              std::vector< double >::const_iterator knot1,
              std::vector< Coord< double > >::const_iterator p0,
              std::vector< Coord< double > >::const_iterator p1,
              std::vector< double >::const_iterator w0,
              std::vector< double >::const_iterator w1,
              size_t p )
    : BSplineCurveBase( knot0, knot1, p0, p1, p )
    { setWeight( w0, w1 ); }

  /// 配列全体を指定して構築
  ///
  /// knot : ノット列
  /// pnt : 点列
  /// w : 重み値の列
  /// p : 次数
  NurbsCurve( const std::vector< double >& knot,
              const std::vector< Coord< double > >& pnt,
              const std::vector< double >& w,
              size_t p )
    : BSplineCurveBase( knot, pnt, p )
    { setWeight( w.begin(), w.end() ); }
};

/*
  NurbsCurve::setWeight : 重み値の取り込み

  w0,w1 : 重み値の範囲 (w1 は末尾の次を表すため含まれない)
*/
void NurbsCurve::setWeight( vector< double >::const_iterator w0,
                            vector< double >::const_iterator w1 )
{
  w_.assign( size(), 1 ); // 重み値のデフォルトは 1 とする

  size_t sz = std::distance( w0, w1 );
  if ( sz <= size() )
    std::copy( w0, w1, w_.begin() );          // データは範囲内
  else
    std::copy( w0, w0 + size(), w_.begin() ); // データ数が超過している場合、末尾は捨てる
}

/*
  NurbsCurve::calc : 曲線の X/Y 座標値を計算する

  pnt : 点列の X/Y 座標
  t : 計算対象の独立変数

  返り値 : 計算結果
*/
double NurbsCurve::calc( const vector< double >& pnt, double t )
{
  size_t p = order(); // 次数

  vector< double > bs( size() ); // B-Spline 基底関数の計算結果(重み値を積算)
  double wSum = 0; // bs の合計
  for ( size_t j = 0 ; j < size() ; ++j ) {
    bs[j] = w_[j] * GCurve::CalcBSplineBasisFunc( knot_, j, p, t );
    wSum += bs[j];
  }

  double ans = 0; // 計算結果
  for ( size_t j = 0 ; j < size() ; ++j ) {
    double d = pnt[j] * bs[j] / wSum;
    if ( ! isnan( d ) )
      ans += d;
  }

  return( ans );
}

NURBS 曲線を描画するためのクラス NurbsCurve も、BSplineCurveBase からの派生になっています。BSplineCurve などとの違いは、構築時に重み値を受け取る箇所と、座標値の計算式のみです。重み値はメンバ変数 w_ に格納しますが、デフォルトは 1 として、引数の範囲が点列の数より小さければ末尾側をデフォルトのままとし、大きければ引数の配列の末尾側を無視するような処理になっています。

NURBS を利用することで、例えば円を簡単に描画することができます。円の方程式は

なので、0 ≤ θ ≤ π / 2 の範囲のみを考えると

より

となります。一方、点列を P₀ = ( 1, 0 ), P₁ = ( 1, 1 ), P₂ = ( 0, 1 ) の 3 点とした時、二次のバーンスタイン基底関数を利用して NURBS を表現すると

となり、w₀ = w₁ = 1, w₂ = 2 とすれば円の方程式と一致することがわかります。

円全体を描画する場合はもう少し複雑で、ノット列を { 0, 0, 0, 1/4, 1/4, 1/2, 1/2, 3/4, 3/4, 1, 1, 1 } とし、重み値を { 1, 1/√2, 1, 1/√2, 1, 1/√2, 1, 1/√2, 1 }、点列を ( 0, 1 ), ( 1, 1 ), ( 1, 0 ), ( 1, -1 ), ( 0, -1 ), ( -1, -1 ), ( -1, 0 ), ( -1, 1 ), ( 0, 1 ) とすることで描画できます。下図は、この条件で B-スプラインと NURBS の両方を使って曲線を描画した結果を示しています。

図 3-1. NURBS による円の描画

B-スプライン	NURBS

NURBS で描画した側は、四隅にある各点に対してより離れた位置に曲線が移動している様子が確認できます。この点に対する重み値が 1 より小さいため、点の位置による影響が小さくなっていることがこの結果からもわかります。重みを大きくすれば影響が大きくなるので、曲線は四角形に近い形となるでしょう。

今回は、B-スプライン曲線の描画方法を中心に紹介しました。ベジェ曲線は B-スプライン曲線の中に含まれることから、さらに多様な形状も表現することが可能となり、ベジェ曲線を使ったスプライン曲線と比べても、より滑らかな曲線を描画することが可能となります。

前に戻る

タイトルに戻る