Bezier曲線による関数の近似 - その5 [レンズ形状をBezierで描く]
Bezier曲線を関数近似に使いたいと言う話の続き。2次元の曲線の話をここでひと段落ということで、実装して精度を確認．．．．
直接解くなら4次方程式を解くことになって、しかしそれは結構大きな式になるので、ヤコビアンを直接書き下してニュートン法で数値的に解けば速い。
今回もまずはMathematicaで書いた。MathematicaのFindRootはかってにヤコビアンを計算してくれるので、面倒が少ない。
3次は誤差0でフィットできるのは当たり前なので、それ以上の4、7、10、12、15次の関数を近似した。実線がベキ関数、点線がBezier曲線による近似。あたりまえだけど結構違っている。
$y$方向の差をプロットするとこうなる。
当然ながら次数が上がるほど近似誤差が増える。差の最大の$x$の位置が次数にしたがって右に移動するのは、$x$が大きくなるほど急峻になって、右下かどが鋭くなるのにBezierは3次なので急には曲がれないからである。
こんなのを考える。 \begin{equation} f(x) = \cos \frac{\pi(6+n)x}{6} \end{equation} これを$n=0\dots$として、半波長から$\pi/6$ずつ伸ばしていく。
実線がもとの余弦関数で、点線がBezier近似の曲線。
$y$方向の差をプロットしたのがこれ。
$(0,\pi)$区間では余弦関数と0.2%ほどの誤差がある。$(\pi,2\pi)$はこの符号が逆なだけなので、繋いでいけば関数全体が近似できる。もちろん、端点での一致や、その微係数の一致を無視すればもっと誤差を減らすことができる。誤差最小だけを狙うのならそっちの方がいい。
しかし、Bezier曲線とは関係ないけど、正弦余弦関数は3次のベキでここまで近似できるということで、これは僕にはちょっと意外だった。
区間が$(0,\pi)$を超えると急速に精度が悪くなる。この悪くなり方はほんとに急速でなんでそんなに違うのかピンときにくい。これはこれで面白い。
4分の1円の場合に、微係数を合わせる方法では、たまたまちょうどキリのいい数字になった（$\sqrt{7}$がキリがいいとしてだけど）。
ここでは中途半端な円弧を考えてみる。
非球面レンズを使ってる光学屋なら、つい円弧の式を
\begin{equation}
y=\frac{x^2}{r \left(1+\sqrt{\displaystyle 1-\frac{x^2}{r^2}}\right)}
\end{equation}
と書いてしまう。この式は普通の人は「これが円弧の式？」というかもしれないけど、これは非常に便利な形。これをそのまま使おう。
これを$x=0 \dots x_0$までの円弧だとして、微分を計算すると \begin{align} \frac{dy}{dx} &= \frac{x}{\displaystyle r\sqrt{1-\left(\frac{x}{r}\right)^2}} \\ \frac{d^2y}{dx^2} &= \frac{1}{\displaystyle r\left(1-\left(\frac{x}{r}\right) \right)^{3/2}} \end{align} となるので、端点$\boldsymbol{P}_s=(0,0)$と$\boldsymbol{P}_e=(x_0,x_0^2/(r(1+\sqrt{1-(x_0/r)^2}))$とでは \begin{align} {}_1C_s &= 0 \nonumber \\ {}_1C_e &= \frac{x_0}{\displaystyle r\sqrt{1-\left(\frac{x_0}{r}\right)^2}} \nonumber \\ {}_2C_s &= \frac{1}{r} \nonumber \\ {}_2C_e &= \frac{1}{\displaystyle r\left(1-\left(\frac{x_0}{r}\right) \right)^{3/2}} \end{align} として式-\ref{c2eq}を一致させるように中間の制御点を決めれば、円弧をBezier曲線で近似できたことになる。
前節のやりかたで4分の1円より小さな円弧をBezier曲線で近似して、その真円からの最大誤差を円弧の内角に対してプロットしたのが下である。
図の赤点が4分の1円の場合で、当然ながら円弧を狭くしていくと誤差は急速に小さくなる。内角が45°、つまり8分の1円の場合は4分の1円に対して100分の1ぐらいに、30°の12分の1円にすると数千分の1にできる。
たとえばIllustratorで円を描くと、4分の1円を4つ、つまり固定した制御点が4つで描くことになるけど、その場合半径1メートルの円を描くと45°方向への半径は2mm短かくなる。
それをこのやり方に従って固定した制御点を8つにすると半径1メートルに対して2$\mu$mほど短いだけになる、ということである。
こう言う誤差を減らすためにBezier関数を拡張して3次以上を使ったり、ベキだけではなく有理関数まで広げた有理Bezierというのがある。でも急速に扱いが難しくなって、はっきり言って僕にはメリットを感じられない。次数を上げれば表現力があがるのは当然だけどそれで得られるメリットが、ここで書いたような0.001以下の誤差を問題にする場合にしか有効ではない。また有理関数は飛躍的な表現力の拡大が得られるけど、数学的な扱いはベキだけの場合と全く違った注意力が必要になる。つまり分母が0になる場合をいつも気をつけていないといけないし、微分や積分も急に難しくなる。
僕としては有理関数にして誤差なしで表現できたからと言って、そりゃそうだろ、としか思えなくて、それよりも3次Bezier曲線を使いこなす方向が性に合ってる、と言う感じ。
レンズは光軸$z$軸のまわりに回転対称なので、体積は \begin{equation} V = \int \pi y^2 \frac{dz}{dy} dy \end{equation} として式-\ref{aspherical}を微分して代入して積分すればいい。
ところが非球面を表す式が一般的に \begin{equation*} {\rm asp}(y) = \frac{\rho y^2}{\sqrt{1-(k+1) \rho ^2 y^2}+1}+a y^4+b y^6+c y^8+d y^{10} \cdots \end{equation*} みたいな格好をしている。数値積分すると不安定になることがたまにあって、たとえば$k$の値のちょっとした違いで、ほとんど形は変わらないのに数値があばれたりする。
なぜ数値的に不安定になるかと言うとベキの項の$y$についての次数が高すぎて、IEEE754のDouble幅（Swiftの癖がついてつい頭を大文字にしてしまう）では項の足し算で桁落ちしてしまうことがあるのが原因である（もう少し正確に言えば、第1項とそれ以降が直交していないので、$k$の値がちょっと変わるとベキの係数がゴロっと変わってしまって計算結果の有効数字が安定しないせい）。しかしもともと光学的にまともな非球面であれば曲率の変化も滑らかで、凸凹はないので本来数値的には安定なはずである。ようするに単にベキ級数を使った表現の問題である。
もっと安定に、サクっと数値積分したかった。Bezier曲線で近似してから積分すると次数が低くて簡単で、しかも体積計算するには十分な精度が出せる。
今回は結構正解だった。満足。
これでひと段落だけど、あとこれを3次元空間内のBezier曲面で考えたい。最終的にはこれを極めたいんだけど、これがまた、急に難しくなるんだわ。ほんとに、もう。
8 テスト
実際にちゃんと動作するかどうか確かめる。 もう一度たどり着いた式をちょっとみやすく書くと \begin{align} \frac{2 Q_s^2 \left(Q_e \left((x_s-y_s+y_e)- {}_1C_s x_e\right) +({}_1C_s-{}_1C_e)\ell_2 \right)}{3 Q_e\ell_1^2} - {}_2C_s &= 0 \nonumber \\ \frac{2 Q_e^2 \left(-Q_s \left( (x_s-y_s+y_e)- {}_1C_e x_e\right)+ ({}_1C_s-{}_1C_e)\ell_1\right)}{3 Q_s\ell_2^2} - {}_2C_e &= 0 \nonumber \end{align} ただし、 \begin{align} Q_s &= \sqrt{{}_1C_s^2+1} \nonumber \\ Q_e &= \sqrt{{}_1C_e^2+1} \nonumber \end{align} で、端点の相対的な位置$(x_s,y_s)$と$(x_e,y_e)$、端点での1階の微係数${}_1C_s$と${}_1C_e$を定数として含んだ式を、２階の微係数${}_2C_s$と${}_2C_e$との差が最小になるように、端点から中間の制御点までの距離$\ell_s$と$\ell_e$を決めればいい、というのがこないだの結論。直接解くなら4次方程式を解くことになって、しかしそれは結構大きな式になるので、ヤコビアンを直接書き下してニュートン法で数値的に解けば速い。
今回もまずはMathematicaで書いた。MathematicaのFindRootはかってにヤコビアンを計算してくれるので、面倒が少ない。
8.1 ベキ関数の近似
まず意地悪問題。本来表現できない4次以上のベキ関数を近似してみる。3次は誤差0でフィットできるのは当たり前なので、それ以上の4、7、10、12、15次の関数を近似した。実線がベキ関数、点線がBezier曲線による近似。あたりまえだけど結構違っている。
8.2 余弦関数の近似
こんどは$\cos$関数を近似してみる。こんなのを考える。 \begin{equation} f(x) = \cos \frac{\pi(6+n)x}{6} \end{equation} これを$n=0\dots$として、半波長から$\pi/6$ずつ伸ばしていく。
$y$方向の差をプロットしたのがこれ。
しかし、Bezier曲線とは関係ないけど、正弦余弦関数は3次のベキでここまで近似できるということで、これは僕にはちょっと意外だった。
区間が$(0,\pi)$を超えると急速に精度が悪くなる。この悪くなり方はほんとに急速でなんでそんなに違うのかピンときにくい。これはこれで面白い。
9 円弧
9.1 Bezier「円」の分割の問題点
「円」をBezierで近似して、それとは別にBezier曲線を分割することを考えた。近似した時点でBezier曲線で「円」ではないんだけど、たとえば$t=1/2$で分割したとき、新しい端点はもとのBezier曲線の誤差の分だけ「円」の上から外れてしまう。大した量ではないので表示だけなら問題になることはない。しかし、半円単独であつかうことはまずなく、たいていは他の部分をBezier近似したものと接続されるはずである。そのとき、わずかなズレなので、どちらかに寄せておけば、つまり一方の端点をもう一方と一致するように無理やり移動させればいいだけである。しかし放っておくと、例えば数値積分なんかで、エラーを起こす可能性がある。4分の1円の場合に、微係数を合わせる方法では、たまたまちょうどキリのいい数字になった（$\sqrt{7}$がキリがいいとしてだけど）。
ここでは中途半端な円弧を考えてみる。
9.2 4分の1円以下の円弧の場合
一般の場合に微係数に制限をつけたので、4分の1円のときと向きをちょっと変える。これを$x=0 \dots x_0$までの円弧だとして、微分を計算すると \begin{align} \frac{dy}{dx} &= \frac{x}{\displaystyle r\sqrt{1-\left(\frac{x}{r}\right)^2}} \\ \frac{d^2y}{dx^2} &= \frac{1}{\displaystyle r\left(1-\left(\frac{x}{r}\right) \right)^{3/2}} \end{align} となるので、端点$\boldsymbol{P}_s=(0,0)$と$\boldsymbol{P}_e=(x_0,x_0^2/(r(1+\sqrt{1-(x_0/r)^2}))$とでは \begin{align} {}_1C_s &= 0 \nonumber \\ {}_1C_e &= \frac{x_0}{\displaystyle r\sqrt{1-\left(\frac{x_0}{r}\right)^2}} \nonumber \\ {}_2C_s &= \frac{1}{r} \nonumber \\ {}_2C_e &= \frac{1}{\displaystyle r\left(1-\left(\frac{x_0}{r}\right) \right)^{3/2}} \end{align} として式-\ref{c2eq}を一致させるように中間の制御点を決めれば、円弧をBezier曲線で近似できたことになる。
9.3 円弧からの誤差
4分の1円のときと同じように真円からの誤差を評価してみる。4分の1円では$t=0.5$のとき半径の誤差として-0.4%ぐらいあった。Bezier曲線の分割では8分の1円まで$t>0.5$はこの最大誤差は変わらない。前節のやりかたで4分の1円より小さな円弧をBezier曲線で近似して、その真円からの最大誤差を円弧の内角に対してプロットしたのが下である。
たとえばIllustratorで円を描くと、4分の1円を4つ、つまり固定した制御点が4つで描くことになるけど、その場合半径1メートルの円を描くと45°方向への半径は2mm短かくなる。
それをこのやり方に従って固定した制御点を8つにすると半径1メートルに対して2$\mu$mほど短いだけになる、ということである。
こう言う誤差を減らすためにBezier関数を拡張して3次以上を使ったり、ベキだけではなく有理関数まで広げた有理Bezierというのがある。でも急速に扱いが難しくなって、はっきり言って僕にはメリットを感じられない。次数を上げれば表現力があがるのは当然だけどそれで得られるメリットが、ここで書いたような0.001以下の誤差を問題にする場合にしか有効ではない。また有理関数は飛躍的な表現力の拡大が得られるけど、数学的な扱いはベキだけの場合と全く違った注意力が必要になる。つまり分母が0になる場合をいつも気をつけていないといけないし、微分や積分も急に難しくなる。
僕としては有理関数にして誤差なしで表現できたからと言って、そりゃそうだろ、としか思えなくて、それよりも3次Bezier曲線を使いこなす方向が性に合ってる、と言う感じ。
10 そもそも何をやりたかったか
なんで突然思いついたようにBezier曲線で近似しようと思ったかというと、実は非球面レンズの体積を計算しようとした。レンズの面をふつう \begin{equation} z = {\rm asp}(y) \label{aspherical} \end{equation} と表す。$y$は光軸からの高さで、$z$は光軸方向にどのくらいずれてるかを表していて、この$z$の値のことをレンズ屋はサグ(sag）と呼ぶ。レンズは光軸$z$軸のまわりに回転対称なので、体積は \begin{equation} V = \int \pi y^2 \frac{dz}{dy} dy \end{equation} として式-\ref{aspherical}を微分して代入して積分すればいい。
ところが非球面を表す式が一般的に \begin{equation*} {\rm asp}(y) = \frac{\rho y^2}{\sqrt{1-(k+1) \rho ^2 y^2}+1}+a y^4+b y^6+c y^8+d y^{10} \cdots \end{equation*} みたいな格好をしている。数値積分すると不安定になることがたまにあって、たとえば$k$の値のちょっとした違いで、ほとんど形は変わらないのに数値があばれたりする。
なぜ数値的に不安定になるかと言うとベキの項の$y$についての次数が高すぎて、IEEE754のDouble幅（Swiftの癖がついてつい頭を大文字にしてしまう）では項の足し算で桁落ちしてしまうことがあるのが原因である（もう少し正確に言えば、第1項とそれ以降が直交していないので、$k$の値がちょっと変わるとベキの係数がゴロっと変わってしまって計算結果の有効数字が安定しないせい）。しかしもともと光学的にまともな非球面であれば曲率の変化も滑らかで、凸凹はないので本来数値的には安定なはずである。ようするに単にベキ級数を使った表現の問題である。
もっと安定に、サクっと数値積分したかった。Bezier曲線で近似してから積分すると次数が低くて簡単で、しかも体積計算するには十分な精度が出せる。
今回は結構正解だった。満足。
これでひと段落だけど、あとこれを3次元空間内のBezier曲面で考えたい。最終的にはこれを極めたいんだけど、これがまた、急に難しくなるんだわ。ほんとに、もう。
2021-02-07 21:10
nice!(0)
コメント(0)
コメント 0