Bezier曲線による関数の近似 - その3

Bezier曲線に対して、スプライン補間のように微係数の値を与えて、しかし補間ではなく関数近似に使おうという試みの続き。特別な場合として4分の1円をやった。今日はもうすこし一般的な場合について。

6  Bezier一般の近似

一般の場合、といっても最初に書いた制限、つまり$y=f(x)$と$x$の陽な形に書けるという制限の上での話。

まず、中間の制御点$(x_1,y_1)$、$(x_2,y_2)$を両端点からの相対位置で表すことにする。 つまり \begin{align} x_{1s} &= x_1 - x_0 \nonumber \\ y_{1s} &= y_1 - y_0 \nonumber \\ x_{2e} &= x_2 - x_3 \nonumber \\ y_{2e} &= y_2 - y_3 \nonumber \end{align} とする。

両端点での1階と2階の$x$座標に関する微係数 \begin{align} \left. \frac{dy}{dx} \right|_{t=0} &= {}_1C_s \nonumber \\ \left. \frac{dy}{dx} \right|_{t=1} &= {}_1C_e \nonumber \\ \left. \frac{d^2y}{dx^2} \right|_{t=0} &= {}_2C_s \nonumber \\ \left. \frac{d^2y}{dx^2} \right|_{t=1} &= {}_2C_e \nonumber \end{align} が与えられているとする。ただしどれも0でもなければ$\infty$でもないとする。これから$x_1$、$x_2$、$y_1$、$y_2$を決める。

端点の微係数の式から \begin{align} \frac{y_{1s}}{x_{1s}} &= {}_1C_s \\ \frac{y_{2e}}{x_{2e}} &= {}_1C_e \end{align} となるので、端点から見た中間の制御点の方向は決まる。これはIllustratorなどを使っていると感覚的によくわかる。

方向は決まるので距離を \begin{align} \ell_1 &= \sqrt{x_{1s}^2+y_{1s}^2} \\ \ell_2 &= \sqrt{x_{2e}^2+y_{2e}^2} \label{lengths} \end{align} とする。こうすると \begin{align} \begin{array}{rl} x_{1s} &= \pm \frac{\displaystyle \ell_1}{\displaystyle \sqrt{1+{}_1C_s^2}} \\ &\cdots \end{array} \label{lengexpr} \end{align} などと書けるが、平方根を開くので復号がついてしまう。

しかし最初にあげた$x$に関して陽に書ける、という制限から \begin{align} x_0 &\le x_1 \nonumber \\ x_2 &\le x_3 \nonumber \end{align} つまり \begin{equation} x_{2e} \le 0 \le x_{1s} \label{explicit} \end{equation} となるので、複号は微係数の符号で決まって \begin{align} x_{1s} &= \frac{\displaystyle \ell_1}{\displaystyle \sqrt{1+{}_1C_s^2}} \nonumber \\ y_{1s} &= \frac{\displaystyle {}_1C_s \ell_1}{\displaystyle \sqrt{1+{}_1C_s^2}} \nonumber \\ x_{2e} &= -\frac{\displaystyle \ell_2}{\displaystyle \sqrt{1+{}_1C_e^2}} \nonumber \\ y_{2s} &= -\frac{\displaystyle {}_1C_e \ell_2}{\displaystyle \sqrt{1+{}_1C_e^2}} \label{lengexpr2} \end{align} とすることができる。

ここのところを工夫して、例えば$dx/dt$と$dy/dt$の符号を$dy/dx$と一緒に与えることにすれば（つまり始点からどの方向に出るかを指定すれば）陽関数であるという制限ははずすことができる。けど今回はもういい。

ただし、一つ問題があって、式-\ref{explicit}は陽に書ける条件だけど、 \begin{equation*} \begin{array}{rl} x_2 &< x_1 \\ \mbox{や、あるいは}& \\ x_2 &\le x_0 \\ \mbox{あるいは}& \\ x_3 &\le x_1 \end{array} \end{equation*} などとなる可能性を排除できない。この場合でも成り立つようにできるかどうかよくわからない。また、微係数を与えるだけなので、上のような条件になっていないかはわからない。どうなるかはやりながら考えることにする。

ということで、式-\ref{lengths}の$\ell_1$と$\ell_2$を2階の微係数から決めることを考える。

2階の微係数は第1回の式-9などだったので \begin{align} 3 {}_2C_s x_{1s}^3+2 \left((x_{2e} + x_{es}) y_{1s} - x_{1s}(y_{2e}+y_{es}+x_0-y_0) \right) &= 0 \\ 3{}_2C_e x_{2e}^3 + 2 \left((x_{1s} - x_{es})y_{2e} - x_{2e}(y_{1s}-y_{es}-x_0+y_0) \right) &= 0 \end{align} ただし、$x_{es}$と$y_{es}$は終点を始点からみた相対座標である。そっくりだけど微妙に符号なんかが違う。あってるかな。

式-\ref{lengexpr}をこの二つの式に代入して$\ell_1$と$\ell_2$に関して解けばいいのだけど、これは巨大な式になる。もとの式に3乗が含まれていて、その中に平方根が含まれるので大きくなるのは当たり前で、前回は少しでも簡単にするために変形した結果4次方程式に落とし込んだ。

今回はそう言うやり方ではなく、少し違うことをする。2次の微係数に関して整理して、その式が2次の微係数に収束するように最適化する。$\ell_1$と$\ell_2$についての2次元の探索になる。

${}_2C_s$と${}_2C_e$を解くと \begin{align} {}_2C_s &= \frac{2 Q_s^2 \left(Q_e \left((x_0-y_0+y_{es})- {}_1C_s x_{es}\right) +({}_1C_s-{}_1C_e)\ell_2 \right)}{3 Q_e\ell_1^2} \nonumber \\ {}_2C_e &= \frac{2 Q_e^2 \left(-Q_s \left( (x_0-y_0+y_{es})- {}_1C_e x_{es}\right)+ ({}_1C_s-{}_1C_e)\ell_1\right)}{3 Q_s\ell_2^2} \end{align} ただし \begin{align} Q_s &= \sqrt{{}_1C_s^2+1} \nonumber \\ Q_e &= \sqrt{{}_1C_e^2+1} \nonumber \end{align} となる。ほんまかいな。

このそれぞれ右辺を${}_2\hat{C}_s(\ell_1,\ell_2)$と${}_2\hat{C}_e(\ell_1,\ell_2)$と書くと、与えられた${}_2C_s$と${}_2C_e$との差をとってその差を最小にすることを考える。つまり \begin{align} E_s(\ell_1,\ell_2) &= {}_2C_s - {}_2\hat{C}_s(\ell_1,\ell_2) \\ E_e(\ell_1,\ell_2) &= {}_2C_e - {}_2\hat{C}_e(\ell_1,\ell_2) \end{align} として \begin{equation} \boldsymbol{E}\equiv \big(E_s(\ell_1,\ell_2),E_e(\ell_1,\ell_2) \big) = \boldsymbol{0} \label{vectorequation} \end{equation} の解を探索すればいい。一見大変そうだけど、$\ell_1,\ell_2$以外は定数なので \begin{align} E_s(\ell_1,\ell_2) &= {}_2C_s - \frac{p_s +q_s \ell_2}{r_s \ell_1^2} =0 \\ E_e(\ell_1,\ell_2) &= {}_2C_e - \frac{p_e +q_e \ell_1}{r_e \ell_2^2} =0 \end{align} のような形を$\ell_1,\ell_2$について解けばいい。ここで$p_s,q_s,r_s,p_e,q_e,r_e$は与えられた微係数から決まる定数である。

これは2変数の2次連立方程式だけど、変数がサイクリックに入っているのでひとつの変数に対して（3次の項のない）4次方程式を解くことになる。4次方程式はいちおう解の公式があるので力づくで解くことはできる（Mathematicaは解の公式を知っていて、Solveで解くと瞬時に返る）。

しかしご存知のように解の公式はすごく大きな式になるのと、解は4つあって、もし実数解があるとすれば4つとも実数か、あるいはふたつが実数で残りは共役なふたつの複素数のはずである（2重解3重解をそれぞれふたつみっつと数えて）が、計算を普通のIEEE754で計算したとすると、丸め誤差や桁落ちで判断の難しい解が得られる場合が十分ある。

つまり、公式があるとは言え、必要な解を間違いなく選ぶのはあんがい難しい。しかもこれでは前回1変数の4次方程式に落とし込んだのと（式は違うけど）まったく同じでつまらない。

直接解かないとすると、式-\ref{vectorequation}を2次元のニュートン法で逐次的に解く方法が考えられる。ヤコビアンも簡単に計算できるので収束は速い。また、あるいは \begin{equation} E = \boldsymbol{E} \cdot \boldsymbol{E} \end{equation} を最小化すればいい。$\cdot$は内積である。これは次数はあがるけど、1次元になるので探索アルゴリズムから実装するにはこっちの方が楽になる。Mathematicaを使えばFindMinimumを使うかFindRootを使うか、の違いでしかなくなるけど。

このどちらにしてもヤコビアンの形を見ればわかるけど、ふたつの変数のそれぞれにローカルミニマムのできにくい形をしている。前回の4次方程式を逐次近似で解くよりずっと速く、しかもずっと安定のはずである。これが今回言いたかったこと。

これを確認するために、続く。