NURBSの微分で躓いた話 - エフアンダーバー

気まぐれで実装したNURBSの微分で何時間も悩んだので備忘録も兼ねて記事に。大概のバグがそうであるように、オチは結構しょうもないのですが・・・

記事中でところどころ歯切れが悪いのは、専門外の内容をネットで拾い集めた知識で書いているからです。間違っている可能性が多分にあるので、鵜呑みにはしないように。

NURBS

とりあえずはNURBSとは何かという話から。

NURBSは Non-Uniform Rational Basis Spline の略で、ざっくり言うと自由度の高い曲線です。複数の制御点といくつかのパラメータを与えることで制御点に沿った形状の曲線を得ることができます。 Non-Uniformは非一様の意味でノットの間隔が一定でないことを表しており、このため各点が影響を与える範囲の調節が可能です。 Rationalは有理の意味ですが、より単純には各点にウェイトを設定でき、各点の影響度を調節できることを表すそうです。 Basisは基底、Splineは複数の部分的な曲線を滑らかにつないでできた曲線を表しますが、単純にB-Spline曲線という曲線の拡張の意と捉えるのがいいかと思います。

ちなみに Uniform Non-Rational Basis Spline （つまり、ノットの間隔とウェイトが一定のNURBS）を考えると、かの有名なベジェ曲線になるそうです。各曲線の関係としてはベジェ曲線の拡張がB-Spline曲線で、さらにその拡張がNURBSと言えそうです。

定義

n次のNURBSの定義は以下の通り。

$\displaystyle \begin{align*} m &: \mbox{制御点の個数} \\ P_{i} &: \mbox{i番目の制御点の座標} \\ w_{i} &: \mbox{i番目の制御点のウェイト} \\ t_{i} &: \mbox{i番目のノット} \end{align*}$

として、

$\displaystyle C(t) = \sum_{i=1}^{m} R_{i,n}(t)P_{i} = \cfrac{\sum_{i=1}^{m} w_{i}N_{i,n}(t)P_{i}}{\sum_{j=1}^{m} w_{j}N_{j,n}(t)} \qquad for \quad t \in [t_{n + 1}, t_{m + 1}] \\ R_{i,k}(t) = \cfrac{w_{i}N_{i,k}}{\sum_{j=1}^{m} w_{j}N_{j,n}(t)}$

ここで $N_{i,k}(t)$ はB-Spline基底関数で、

$\displaystyle N_{i,0}(t) = \left\{ \begin{array}{cl} 1 & if \quad x_{i} \le t \lt x_{i+1} \\ 0 & otherwise \end{array} \right. \\ N_{i,k}(t) = \cfrac{t - t_{i}}{t_{i+k} - t_{i}}N_{i,k-1}(t) + \cfrac{t_{i+k+1} - t}{t_{i+k+1} - t_{i+1}}N_{i+1,k-1}(t)$

ただし、右辺の分母が0になる場合にはその項を0として扱います。

De Boorのアルゴリズム (B-Spline)

NURBSの微分の話に行く前に少し寄り道。

B-Spline曲線上の点を求めるアルゴリズムとしてDe Boorのアルゴリズムというものがあります。このアルゴリズムは少しの変更でNURBSに対しても適用できます。 B-Spline曲線、NURBSはともに定義をそのままプログラムとして書くと、再帰やら誤差やらで面倒なことになります。その辺をうまいことどうにかしてくれるのがDe Boorのアルゴリズムになります。

ちなみにB-Spline曲線の定義は以下の通りで、NURBSのウェイトを一定にしたものと等しくなります。

$\displaystyle C(t) = \sum_{i=1}^{m} N_{i,n}(t)P_{i} \qquad for \quad t \in [t_{n + 1}, t_{m + 1}] \\$

De Boorのアルゴリズムは以下の通りです。

${t \in [x_{l}, x_{l+1})}$ のとき、

$\displaystyle P_{i,0} = P_{i} \qquad for \quad i = l-n, \ldots, l \\ P_{i,k} = (1 - \alpha_{i,k})P_{i-1,k-1} + \alpha_{i,k}P_{i,k-1} \qquad for \quad k = 1, \ldots, n; \; i = l-n+k, \ldots, l \\ with \quad \alpha_{i,k} = \cfrac{t - t_{i}}{t_{i+n+1-k}-t_{i}}$

と $P_{i,k}$ を定義して、

$\displaystyle C(t) = P_{l,n}$

注目すべきは $P_{l-n}$ から $P_{l}$ までの点しか参照していないということ。定義をよく観察するとわかりますが、それ以外の点は対応する $N_{i,k}$ の値が0になるため影響しません。またDe Boorのアルゴリズムは $k = 0$ から計算していくことで最大長 $n + 1$ のバッファで計算できます。

De Boorのアルゴリズム (NURBS)

NURBSにDe Boorのアルゴリズムを適用するために、NURBSとB-Spline曲線の式を見比べてみます*1。

$\displaystyle C_{NURBS}(t) = \cfrac{\sum_{i=1}^{m} w_{i}P_{i}N_{i,n}(t)}{\sum_{i=1}^{m} w_{i}N_{i,n}(t)} \\ C_{B-Spline}(t) = \sum_{i=1}^{m} P_{i}N_{i,n}(t)$

$C_{B-Spline}(t)$ が求まるなら $C_{NURBS}(t)$ を求められそうな気がする・・・はず。

というわけで、ウェイトおよび各制御点の座標にウェイトを乗算した新たな制御点に対してDe Boorのアルゴリズムを適用します（細かい部分は省略）。

$\displaystyle w_{i,k} = (1 - \alpha_{i,k})w_{i-1,k-1} + \alpha_{i,k}w_{i,k-1} \\ w_{i,k}P_{i,k} = (1 - \alpha_{i,k})w_{i-1,k-1}P_{i-1,k-1} + \alpha_{i,k}w_{i,k-1}P_{i,k-1} \\ with \quad \alpha_{i,k} = \cfrac{t - t_{i}}{t_{i+n+1-k}-t_{i}}$

これらを用いて、

$\displaystyle C_{NURBS}(t) = \cfrac{\sum_{i=1}^{m} w_{i}P_{i}N_{i,n}(t)}{\sum_{i=1}^{m} w_{i}N_{i,n}(t)} = \cfrac{w_{l,n}P_{l,n}}{w_{l,n}} = P_{l,n}$

これでDe BoorのアルゴリズムによってNURBS上の点を計算することができました。

計算結果は $P_{l,n}$ とシンプルになっていますが、実際には $P_{l,n}$ を直接求められるわけではないので、 $w_{l,n}P_{l,n}$ と $w_{l,n}$ をそれぞれ求めて一つ手前の計算をすることになります。このとき、 $P_{i}$ を一次元増やして $\hat{P}_{i} = [ \, P_{i} \; 1 \, ]$ としてからウェイトをかけて演算すると二つの値が同時に計算できます。

NURBSの微分

何故そうなるかといった難しい話は論文等にまかせるとして、B-Spline基底関数の導関数は以下のようになるそうです。

$\displaystyle N'_{i,k}(t) = \cfrac{k}{t_{i+k} - t_{i}}N_{i,k-1}(t) - \cfrac{k}{t_{i+k+1} - t_{i+1}}N_{i+1,k-1}(t)$

NURBSの導関数はこれを用いて、

$\displaystyle C'(t) = \cfrac{\sum_{i=1}^{m}w_{i}P_{i}N'_{i,n}(t)\sum_{i=1}^{m}w_{i}N_{i,n}(t) - \sum_{i=1}^{m}w_{i}P_{i}N_{i,n}(t)\sum_{i=1}^{m}w_{i}N'_{i,n}(t)}{(\sum_{i=1}^{m}w_{i}N_{i,k}(t))^{2}}$

となるらしいのですが、さすがにこれを展開してどうにかするのは厳しい。

そう思ってもう少し調べてみると、De Boorのアルゴリズムを使った式がありました。

$\displaystyle C'(t) = n \cdot \cfrac{w_{l-1,n-1}w_{l,n-1}}{w_{l,n}^{2}} \cdot \cfrac{P_{l,n-1} - P_{l-1,n-1}}{t_{l+1} - t_{l}}$

これでNURBSの微分値が効率的に計算できる・・・はずでした。

何に躓いたか？

今回の記事でもそうでしたが、NURBSの話をしているとどうしてもまずB-Spline曲線の話がでてきます。そして、様々な記号が登場してくることで頭がこんがらがってきます。そうなってくると、やってしまうのが関数の混同です・・・

本記事中にまったく同じ記号で別の関数を表していた部分があったのですが、気づいたでしょうか？
それはB-Spline曲線の計算に使用した $P_{i,k}$ とNURBSの計算に使用した $P_{i,k}$ です。もちろん、最後の式にでてくる $P_{i,k}$ はNURBSのほうのものです。が、二つが別のものなどとは露程にも思わず、プログラムを記述するときにB-Spline曲線のほうの式で計算していました・・・

まさか理論レベルで間違っているとは思わなかったので、添え字がずれていないかや、大きな計算誤差が発生する箇所がないかなどと見当違いなチェックを延々と繰り返すはめになりました。そもそも違いがかなり微妙なので、数値微分で近似値を出していなければ気付きすらしなかったかもしれません。テスト大事。

ちなみにちゃんとNURBSのほうの式を用いて計算した結果、正しい値が算出できました。