Bezier曲線による関数の近似

ずっと以前、レンズの断面形状の表示にBezier曲線を使いたくて、非球面をBezier曲線で近似することを考えた。これはそれなりにうまくいったけど、4次方程式を解く必要があって、Mathematica以外では場合分けなんかがめんどうだった。ちょっと面白いやりかたを思いついたので脇道．．．．

光学の分野では1$\mu$m以下のエラー（誤差）が問題、というかその一桁下でも致命的になることがあるので、近似曲線や曲面を使うことは少ないんだけど、一旦光学を離れると1$\mu$mのエラーは無視できることが多い。特に以前やった表示の場面では、光学の問題は原理的に滑らか（波動方程式で記述される場）なので、近似精度よりも曲線曲面が滑らかで曲率の変化や変曲点の存在などが的確に伝わることの方が重要である。

そしてBezierで表現されているといろいろ便利なことが増えてきた。コンピュータグラフィクスの分野ではBezierでなければいきなりもうビットマップしかない、というようにBezierとそれ以外と言っていいほど多く使われるようになっている。

また、外部のソフトウェアへのインターフェイスとしても使われることが増えたし、さらに例えば面積体積の計算などは光線追跡ほどの精度は必要ないけど、非球面式では積分がめんどくさい、また単なる点列やそのスプライン補間では点数が多いのでもっとめんどくさい、というような場面で使える。

Bezier曲線で的確に近似されていれば、少ない点数で振動のない素直な形が表現できて、式も簡単（3次のベキ）なので積分も実行不能になることはない、という利点がある。

ところが、そもそもBezier曲線は滑らかな任意形状を直感的に表現するためのもので、特に図形的なユーザインターフェイスの上で制御点を直感的に扱ってユーザの思い通りの曲線を作ることができる、というのがもっとも有用な点なのでその逆、ユーザが直接制御点を操作するのではなく、何か他の要因から制御点を決定しようというのはほとんどされていない。

しかし思い返せば、Bezier曲線の定義に使われているBernstein関数$b_{n,i}(t)$というのは、もともと多項式による連続関数の近似に使われたものだった。

$b_{n,i}(t)$は区間$(0,1)$で定義された任意の連続関数$f(t)$を近似できる。つまり \begin{equation} B_n(f:t) = \sum_{i=0}^nf(i/n)b_{n,i}(t) \end{equation} とすると、$B_n(f)(t)$は \begin{equation} \lim_{n\rightarrow \infty} B_n(f:t) = f(t) \end{equation} に収束し、しかも一様収束する。実際に低次の和を計算してみると、速くはないがすなおに収束する（一方向から近づく）ことがわかる。独立かつ稠密だけど直交していない関数列のいいところが現れた振る舞いをする。この特徴はすごく面白くて、数値的なフィッティングの問題にも使えるかと思ったんだけど、難しかった。ポイント数と同じ次数のベキ関数になるので、ちょっとポイント数が多いとIEEE754のDouble精度でも足りなくなってしまう。残念。

Bezier曲線は3次なのでこの一様収束性などの特徴が使えるわけではないが、素直さは役に立ちそうな気がする。そこでもういちど掘り返してもう少し一般的な場合に簡単にできる方法を考えたい。

2  Bezier曲線の定義

Bernstein関数$b_{n,i}(t)$の定義は \begin{equation} b_{n,i}(t) = {\binom ni} \;t^i \, (1-t)^{n-i} \end{equation} で、$ {\binom ni}$は2項係数である。 Bezier曲線では3次まで使うことになっていて \begin{align} \boldsymbol{B}(t) &= \sum_i \boldsymbol{P}_i b_{n-1,i}(t) \\ &= \sum_{i=0}^3 (x_i,y_i) b_{3,i}(t) \end{align} ここで$\boldsymbol{P}_i = (x_i,y_i)$は2次元の点の位置で、制御点と言われる。制御点は3次元の空間座標で使われることもある。

Bezier曲線の性質などは、ここよりも他所様の親切なサイトを確認した方がよろしいです。

3  近似の条件

このBezier曲線を関数の近似に使いたいわけだけど、どういう条件で近似するかというと、誤差最小などではかなり難しくなってしまう。

関数全体というよりはある特定の部分を近似して、それをつないでいければいい。ちょうどスプライン補間をするようなイメージで、しかしスプライン補間に比べてずっと少ない制御点で近似することを考える。というか、そもそもスプライン補間は離散点の「補間」が目的であって近似ではない。スプライン補間のような手続きで滑らかに「近似」する曲線を得たい、というのが僕の目的である。

スプライン補間をするときと同じような考え方で両端点の座標と、そこでの1階と2階の微係数を一致させることを考える。ようするに4つの制御点のうち$\boldsymbol{P}_0$と$\boldsymbol{P}_3$は直接決めて、$\boldsymbol{P}_1$と$\boldsymbol{P}_2$を微係数から決めよう、というものである。

また、簡単のため近似対象の関数について条件をつけておく。

1. 連続である
2. $y=f(x)$と表せる。つまり多価であってもよいが$x$について陽（explicit）に書ける
3. $x$について2階微分可能である、また両端点で2階微分が0ではない

とする。一つ目は当然だけど、二つ目は、例えば半円は対象にするけど閉じた円はしない、ということで、Bezier曲線のAffine不変性のメリットを生かさない、ということになってしまう。また、三つ目の後半は端点で曲率0ではない、ということで、かなり制約をすることになるけど、その場合はより簡単な手法が使えるのでいいだろう。

とりあえず簡単のためこうしておく。

4  座標に関する微分

そのためにBezier関数の座標に関する（媒介変数の$t$ではなくて）微分を計算しておく。

1階微分は \begin{align} \frac{dy}{dx} &= \frac{dy}{dt} / \frac{dx}{dt} \\ &= \frac{y_0 (t-1)^2+y_1 (3t-1)(1-t) + y_2(3t-2)t + y_3 t^2}{x_0 (t-1)^2+x_1 (3t-1)(1-t) + x_2 (3t-2)t + x_3 t^2} \nonumber \end{align} なので、端点での値は \begin{align} \left. \frac{dy}{dx}\right|_{t=0} &= \frac{y_0-y_1}{x_0-x_1} \nonumber \\ \left. \frac{dy}{dx}\right|_{t=1} &= \frac{y_2-y_3}{x_2-x_3} \label{firstdiff} \end{align} となる。当然$x_1=x_0$などだったりすると計算できないことになるけど、それはここでは考えず、そういう場合は別途他の手法を考えることにする。

2階微分は \begin{equation} \frac{d^2y}{dx^2} = \frac{d}{dt}\left(\frac{dy}{dx}\right) / \frac{dx}{dt} \end{equation} と書けるけど、具体的に式を書き下したらかなり大きくなってしまうのでここに書くことはやめる。

しかし端点での値はそれほど大きくならなくて \begin{align} \left. \frac{d^2y}{dx^2}\right|_{t=0} &= \frac{2 \left(x_2 \left(y_0-y_1\right)+x_0 \left(y_1-y_2\right)+x_1 \left(y_2-y_0\right)\right)}{3 \left(x_1-x_0\right){}^3} \nonumber \\ \left. \frac{d^2y}{dx^2}\right|_{t=1} &= \frac{2 \left(x_3 \left(y_2-y_1\right)+x_2 \left(y_1-y_3\right)+x_1 \left(y_3-y_2\right)\right)}{3 \left(x_2-x_3\right){}^3} \label{seconddiff} \end{align} となる。これがわかれば近似はできる。この式-\ref{firstdiff}と式-\ref{seconddiff}を数値的に与えて$(x_1,y_1)$と$(x_2,y_2)$が決定できればいい、ということになる。

しかし、4変数の3次連立方程式なので、力づくで直接解けなくはないが、非常に大きな式になってしまう。前回はいろいろな条件をつけて最終的にひとつの変数に対する4次方程式に帰着させた。

今回は別のアプローチを考える。