現在、BBSに書き込みキーを設けています。「書き込みキー」欄に、「MYOMOTO」をカギカッコは抜いて、 半角、小文字で 打ち込んで書き込みをしてください。 このキーワードは時々変更されますが、その都度こちらにて報告します。
PCを付けて時刻を見ると8時23分、8時にかけたはずのアラーム鳴らなかったな……でもゴミ出しに行かないと……4月だというのにまだ寒いなとぼやきながら、慌てて服を着てゴミ袋を出して戻り、しばらく昨日書いたコードを見直しているとアラームが鳴りだしました。目覚ましに使っているAndroid端末の時刻には8:00の表示が……117を聞いたら確かに8時だとかで、どうやらtime.windows.comにつながらなかったようです。それにしてもなんで1時間進んだりするかなあ？
今回からレイマーチングに戻りました。というわけで今回のネタは……
MandelbulbもMandelboxもここでは公開してませんが前にやったので4Dジュリア集合に挑戦してみました。春らしい桜色のジュリア集合が描画されます。桜の季節だから桜色なのであって他意はないんですよ？
なんとなく始めたのですが、それにしても物凄く栄養価の高い内容でした。具体的に使っているテクニックといいますか、数学的な方法としては
うわぁー、全部書くの大変そうだよ！一つで一回のひとりごち相当の分量が必要そうですね。特にデュアルクォータニオンの扱い方は調べても詳しくは出てこなかったので自前で適当に書いていますからあんまり自信が無いです……週末をずっと頭を抱えて過ごしたものの、出来上がるのはなんか気持ちの悪い図形というのが浮かばれませんが、ここから応用するのはキミだ！
自分では苦労したつもりなので、忘れないようにメモとして残しておきます。今後何かの役には立つかもしれませんが走り書きなので多分に間違いも含まれているかと思いますから、僕以外の方が読まれる際には慎重に読み進めてください。
まずは4Dジュリア集合の成り立ちについて語りましょう。
そもそもなんで三次元でいいのに四次元なのかというと、クォータニオンを使って座標を計算するためです。三元数は作れないので四次元で計算して、レンダリングするときに1つの軸を捨てちゃいます。計算自体は大体2Dのジュリア集合と同じようなノリで行けます。
\[ Q_0 = 任意, Q_{n+1} = Q_n^2 + C \]というクォータニオンの数列\(Q_n\)を考えて、\(\lim_{n \to \infty}|Q_n|\)が\(\infty\)に発散しないような\(Q_0\)の集合が四次元ジュリア集合です。\(C\)もクォータニオンの定数ですからこれも適当な値を入れます。
今回は四次元で計算しますから、\(Q = w + x\vec{i} + y\vec{j} + z\vec{k}\) と書くことにします。大文字と小文字の区別に気を付けて。 \(\vec{i},\vec{j},\vec{k}\)が虚軸の基底ベクトルで、\(w,x,y,z\)がスカラーです。 GLSLのvec4が\((x,y,z,w)\)の順で記述するのでちょっとコードに起こしづらくなりますが、大体の文献で\(w,x,y,z\)の順で書いてあるので不本意ながらここでもそう書きます。
クォータニオンの基底ベクトル同士の乗算についておさらいしておくと
\[ \vec{i} \vec{i} = \vec{j} \vec{j} = \vec{k} \vec{k} = \vec{i} \vec{j} \vec{k} = -1\\ \vec{i} \vec{j} = - \vec{j} \vec{i} = \vec{k}\\ \vec{j} \vec{k} = -\vec{k} \vec{j} = \vec{i}\\ \vec{k} \vec{i} = -\vec{i} \vec{k} = \vec{j} \]乗算すると違う軸の話になる、というのはテンソルの乗法を知らないと混乱しますが、こうなってます。後で取り上げるクォータニオン二重数を自分でいじると結構使う大事な知識です。
まあ二次元が四次元になっただけ、vec2がvec4に変わっただけですからここまでは難しい話ではないです。とりあえず、\(Q_0\)が任意の数からスタートするので、安直にワールド座標\((x,y,z)\)をそのまま\(Q_0.xyz\)に入れ、\(Q_0.w\)には適当な定数を入れて描いてみました。
……なんだかつまらない図形しかできないんです。\(w\)を固定してしまうと必ず回転体になってしまうんですね。なので、\(Q_0\)の\(x\)座標を適当な定数で固定し、ワールド座標\((x,y,z)\)を\(Q_0(y,z,w)\)に入れて描画する事にしました。つまり、四次元ジュリア集合の、ある\(x\)における断面図が表示されるわけです。「一定の\(w\)における断面図は回転体になってしまう」という点は自分で作るとひっかかるポイントその1だと思います！なお、\(C\)の方は固定でOKのようです。今回のデモのようにフレームごとにグリグリ変形させたい場合には\(C\)をフレーム毎に変化させていきます。
逆に\(Q_0(x) += sin(y);\) という感じで\(Q_0\)のx軸もフレーム内で変化させてみると、こんな感じに複雑さが増します。sin関数のせいでいかにも四次元！て感じでゆらめいていますねぇ！こんなのが 「\(Q_{n+1} = Q_n^2 + C\)として\(Q_n\)が発散しない\(Q_0\)の集合」 という一文だけで定義できちゃうのが面白い所ですね。
ちなみに、3D空間上の座標を\(y,z,w\)に入れたという事を忘れていると法線を計算するときにドはまりします。僕はこれで丸一日と半日ムダにしました。
次に、どうやってジュリア集合を描くのか、という話に進みます。このトピックが今回の核心であり、将来に亘って最も役に立ちそうな内容だと思います。
2Dの場合だったら画面上の座標を適当に掛け算足し算して\(Z_0\)にブチ込んでそれだけの話ですから何も悩むことはないでしょう。慣れている方なら点を打つ機能さえ使えれば覚えたての言語でもいきなり書ける程度と思いますが、3D以上の場合はどうしたらいいんでしょうか？ボリュームレンダリングでもします？
さすがに2017年現在、Webブラウザでフラクタル図形のボリュームレンダリングをおっぱじめられたら皆逃げ帰るでしょう。何事にも工夫が必要です。そこで凄い味方が居ます。Hubbard-Douady potentialです。
Hubbard-Douady potentialというのは以下の\(G(\vec{p})\)で与えられる「\(G(\vec{p})\)が0になるような(複素数空間内での)位置ベクトル\(\vec{p}\)の所にマンデルブロ集合やジュリア集合の境界が出来るよ」と教えてくれる頼もしい奴です。
\[ G(\vec{p}) = \lim_{n \to \infty} \frac{1}{2^n} \log |Q_n| \]※\(\vec{p}\)はジュリア集合の場合は初期値\(Q_0\)マンデルブロ集合の場合は\(C\)の事です。
まずマンデルブロ集合でもジュリア集合でも同じので行けちゃうというのが凄すぎますが、どこからこんなひみつ道具が出てきたのかという話はここでは割愛します。
ともあれ、\(G(\vec{p})=0\)になった所に面を張ればいいんです、やったね！
……どこだよ！それは！と思うでしょう。僕もそうでした。そこで、「等値面のレンダリングをする」ためのテクニックが必要になるのです。
そのための方法は、iq大明神の受け売りになるのですが、\(G\)のようにベクトルを入れるとスカラーを返す、いわゆるスカラー場\(f\)について
\[ f(\vec{p}+\vec{\varepsilon}) = 0 \]として、\(\vec{\varepsilon}\)は十分短いベクトルを考えます。\(\vec{p}\)は等値面表面まで"わりと惜しい"位置にあると考えて……つまり、着目している\(\vec{p}\)から\(\vec{\varepsilon}\)だけレンダリングしたい等値面まで離れているのです。
テーラー展開すると
\[ f(\vec{p}+\vec{\varepsilon}) = f(\vec{p}) + \nabla f(\vec{p}) \cdot \vec{\varepsilon} + \cdots \]\(\nabla f(\vec{p})\)というのは\(\mathrm{grad} f(\vec{p})\)です。\(f(\vec{p})\)を各方向に偏微分して得た勾配ベクトルです。解らなかったら「ベクトル解析学 基礎」とかで調べてみてね！
さて、\(\vec{\varepsilon}\)は小さいので\(\nabla^2\) 以降の項はほぼ無視して
\[ 0 = f(\vec{p}+\vec{\varepsilon}) \fallingdotseq f(\vec{p}) + \nabla f(\vec{p}) \cdot \vec{\varepsilon} \]そしたら、内積やら絶対値やらの性質から以下を得ます。
\[ 0 \ge |f(\vec{p})| - |\nabla f(\vec{p}) \cdot \vec{\varepsilon}| \ge |f(\vec{p})| - |\nabla f(\vec{p})||\vec{\varepsilon}| \]こんな事よく考え付くなぁと感心するしかないのですが、さらに変形して以下を得ます(やってみよう！)
\[ |\vec{\varepsilon}| \ge \frac {|f(\vec{p})|} {|\nabla f(\vec{p})|} \]距離を与えてくれるわけでもない変な関数(さっきまで頼もしい凄い奴って言ってたのに……)しか分からなかったのに「少なくともあとこれだけ行ったら等値面です」って情報が分かりました！距離さえわかればレンダリングできるというのはレイマーチャー(造語)なら本能的に悟れますから、一気にレンダリングへの道が開けましたよ！マクガイバーになった気分ですねー
さて、ここでHubbard-Douady potentialに戻りましょう。
\[ G(\vec{p}) = \lim_{n \to \infty} \frac{1}{2^n} \log |Q_n| を微分して\\ |G'(\vec{p})| = \lim_{n \to \infty} \frac{1}{2^n} \frac{|Q'_n|}{|Q_n|} \]よって、\(\vec{p}\)から等値面までの距離を\(d\)とすると
\[ d \ge \frac{G(\vec{p})}{|G'(\vec{p})|} = \lim_{n \to \infty} \frac{|Q_n| \log|Q_n|}{|Q'_n|} \]コードにした場合どうなのか、というのはこのページのソースにGLSLのコード丸ごと入ってますから、その中のscene関数を読んでみてください。この文章もほぼソースコードのコメントからコピペしている事が丸わかりですね！
等値面のレンダリングの式の中にも厄介な奴が居ることはお分かりでしょう。そうです。\(Q'_n\)です。\(Q_n\)は級数を収束/発散させる事で求めるのにどうやって微分するんでしょうか。また、法線を求めるのにも\(Q'_n\)は有用です。
一般にレイマーチング法では法線の計算には距離関数をてきとーに数値微分して求めます。ですよね？ともかくDE(vec3 p)という関数で点\(\vec{p}\)から物体の表面までの距離が拾える関数があるとしたらこんな感じで
vec3 getNormal(vec3 p){
float d = 小さい数;
return normalize(vec3(
DE(p + vec3( d, 0.0, 0.0)) - DE(p + vec3( -d, 0.0, 0.0)),
DE(p + vec3(0.0, d, 0.0)) - DE(p + vec3(0.0, -d, 0.0)),
DE(p + vec3(0.0, 0.0, d)) - DE(p + vec3(0.0, 0.0, -d))
));
}
数式で書くと \(\nabla DE(\vec{p}) = \vec{N}\) (\(\vec{N}\):法線ベクトル) となります。なんだかんだ言って数式は綺麗ですねえ。
こういうのを差分法っていうんですか？忘れましたから「てきとーな数値微分」呼ばわりしています……別に悪い方法でもないですし、見栄えの点で言うと往々にして閾値とか諸々のパラメータの調整の方が計算法そのものより大事になりますが、この方法では法線を求めるために6回追加でレイを飛ばさないとならないので速度の面でもかなり不利になってきます。そんなわけで今回のような等値面モノでなくとも自動微分を使って偏微分方程式の数値解を求めるテクニックは役立つはずです。
今回は使っていませんが光源が少ないシーンならば数値微分するにしても方向微分を使う、というテクニックがあります。これもiq大明神の受け売りですが、
//light : 平行光源ベクトル
float getDiffuse(vec3 p,vec3 light){
float d = 小さい数;
vec3 d_light = d * light;
return normalize(vec3(
DE(p + d_light) - DE(p - d_light));
}
要するに、光線ベクトル沿いに微分しちゃうんですね。これなら2本レイを飛ばすだけでOKだし、わざわざ法線ベクトルを得てから内積を取る、という事もしなくて良いです。うん、まあこれで大体事足りる気もするんですけど……今回はCook-Torranceシェーディングもしますし、折角なので今回は自動微分で法線も正確に得ましょう！
さて、自動微分そのものはC++みたいに演算子の多重定義のできる言語なら本当に「自動」でできてしまって羨ましいのですが、GLSLだと「わりと手動微分」になります。まあ何はともあれ、自動微分のための下準備として、二重数について考えます。
まず、導関数を知りたくても知れないそんな難しいお年頃の関数\(f(x)\)を考えます。そしたら次に、\(x+\epsilon\)として\(x\)に\(\epsilon\)という小さな数をくっつけてみます。とても小さな数なので\(\epsilon ^2 = 0\)として扱ってよいとします。こういう\(1 + 2\epsilon\) , \(x + x'\epsilon\)みたいに\(\epsilon\)の付いた数を二重数(dual number)と言います。\(1 + 2\epsilon\)を\( \langle 1,2 \rangle \)とか\(x + x'\epsilon\)を\( \langle x,x' \rangle \)という感じで表記したりもします。
では、\( f(x) = x^2\) の場合を例に考えてみます。まだ素直だったあの頃、\(f '(x)=2x\)である事は一目瞭然ですが、さっきの\(x+\epsilon\)を使って
\[ f(x+\epsilon) = x^2 + 2x\epsilon + \epsilon^2 = x^2+2x\epsilon \]ふむ、じゃあ、\(g(x) = x^3+2\)ならどうでしょう。
\[ g(x+\epsilon) = x^3 + 3x^2\epsilon + 3x\epsilon^2 + \epsilon^3 + 2 = x^3 + 3x^2\epsilon + 2 \]お分かりいただけたであろうか？なんと、\(\epsilon\)の付いている項が導関数なのである！とはいえこのままではコードから上手い事使えませんから、もうちょっと実戦に即して考えましょう。
僕らが今格闘している関数って漸化式を使って級数の収束先を求める、とかそういう「\(f(x)\)=なんとか」って形で書くことすらまともに出来ない物じゃないですか、それなのに\(f'(x)\)を求めろとか無茶ぶりなわけです。しかし、\(f'(100)\)みたいに特定の\(x\)についての値を求めることすなわち「偏微分方程式の数値解を求める」事は出来るのです。
実際のジュリア集合の場合のコードはソースを見ていただくとして、「\(f(x)\)という関数自体は綺麗に書けなくても特定の\(x\)についての\(f'(x)\)を求める」方法ついて説明します。まず、float func(x) という関数を考えます。func(x)は以下のようなアルゴリズムだとしましょう。
float func(float x)
{
float r = x;
for (int i=0; i<3; i++) {
r = r * 2.0;
r = r + 1.0;
}
return r;
}
xからスタートしてひたすら2を掛けては1を足しまくる関数です。 例では検算できるようにと\(func(x)=(((x*2+1)*2+1)*2+1 = 8x + 5\)ですから手計算でもこの通り行けますけどループの回数が100回とかになったらもう無理でしょう？さて、\(func'(x)\)は求まるのでしょうか？それには、\(\langle r, d\rangle\)という二重数を考えます。
vec2 func_d(float x)
{
float r = x; //初期条件：r = x, d = dx/dx としてスタート
float d = 1.0; //<r=u,d> → <u,du/dx> (= は代入演算子)
for (int i=0; i<3; i++) {
r = r * 2.0; d = d * 2.0; //<r,d> * u = <ru, du>
r = r + 1.0; //<r,d> + u = <r+u, d>
}
return vec2(r,d);
}
入力を浮動小数点数xとし、出力は二重数で返しています。funcをちょっと改造したfunc_d関数に\(x\)を入れるとあら不思議、返り値に(func(x),func'(x))が帰ってきちゃうのです。嘘だと思うなら検算してみよう。
こんな感じで、ある関数を導関数の値まで求まるように改造する手順としては
これだけ！上記のコードのコメントにあるように、二重数に対して実数倍したら、ε部にも実数倍をしています。一方で実数の加算ではε部には何もしません。先ほど「C++だと本当に自動で微分できる」と書いたのは二重数クラスを作って演算子のオーバーロードをすればいいからなのです。とはいえ、それの出来ないGLSLでもさほど難しくない作業で偏微分方程式の数値解もついでに求めることのできるコードに変えられます。
二重数への何らかの演算に対する具体的な操作については英語版Wikipediaのこちらのページが詳しいです。日本語版Wikipediaにも同じところへリンクが張ってあるので諦めて数式くらい英語が混ざってても読みましょう！
さて、以上を踏まえていよいよジュリア集合のための計算を始めてみましょう。
\(G'(\vec{p})\)を求めるにはクォータニオンの二重数を用いなければなりません。とはいってもジュリア集合でやってる計算って\(Q_{n+1} = Q_n ^2 + C\)だけですから、楽勝楽勝、余裕ッスよ！まずは\(Q_n\)を二重数\(\langle Qr_n,Qd_n \rangle\)とおいて
\[ Q_0 = \vec{p}\\ Qr_0 = Q_0\\ Qd_{0(w)} = \vec{1}, Qd_{0(x)} = \vec{i}, Qd_{0(y)} = \vec{j}, Qd_{0(z)} = \vec{k} \]……どうしたんですか、僕の顔に何かついてます？ともかく\(Q_{n+1}\)以降は以下のようになります
\[ Qr_{n+1} = Qr_n * Qr_n + C\\ Qd_{n+1(w)} = Qr_n * Qd_{n(w)} + Qd_{n(w)} * Qr_n,\\ Qd_{n+1(x)} = Qr_n * Qd_{n(x)} + Qd_{n(x)} * Qr_n,\\ Qd_{n+1(y)} = Qr_n * Qd_{n(y)} + Qd_{n(y)} * Qr_n,\\ Qd_{n+1(z)} = Qr_n * Qd_{n(z)} + Qd_{n(z)} * Qr_n,\\ \]\(Qd\)が分裂してしまいました！しかも4つに！ど、どうしたらいいんだ！でも良く見ると4つとも初期値以外は同じ格好をしています。(表記法に自信が無いので正しい表記をご存知の方は教えてください)
これはどういうことかというと、\(Q_n\)はクォータニオンを引数にとるので、いわゆる多変量関数という扱いなのですが、多変量関数の微分にあたっては、各々のパラメータについての偏微分を考慮しないといけないためです。このため、\(Qr_n\)は1つでいいのですが、\(Qd_n\)は\(w,x,y,z\)各方向の単位ベクトルから出発して、各方向の偏微分値を計算していきます。\(\vec{1}\)というのはクォータニオンのw要素の基底ベクトルである実数軸の事です。つまり、\(Qd_{n(x)} = \frac{\partial Q_n}{\partial x}\)であり、\(\vec{p}\)をx方向に微小距離δだけずらしたら、\(Qn\)は\(\delta Qd_{n(x)}\)だけ動きますよ、という事なんです。かくして、\(Qd_n\)は分身の術によって各軸方向を手分けをして旅することになりました。
クォータニオンは掛ける方向によって結果が違うのでそこにも注意しましょう。ともかくこれで\(n\)をどんどん増やして行けば \(\lim_{n \to \infty} Q_n(\vec{p})\)も \(\lim_{n \to \infty} Q'_n(\vec{p})\)も(近似解を)求めることが出来ます。
次は二重数から\(Q_n\)，\(Q'_n\)に戻す方法について……\(Q_n\)の方はすぐ察しが付くと思います、\(Qr_n\)そのままです。が、4つに分かれた\(Qd\)をまとめて\(Q'_n\)を取り出す方法などなど、いいちこをやりながら手計算で書いてますから幾らか間違ってるかもしれません。間違っていたら僕のせいですから僕にツッコミは入れてください！
わざわざ自動微分を使ってまで求めたかったのは\(|G'(\vec{p})| = \lim_{n \to \infty}\frac{1}{2^n}\frac{|Q'_n|}{|Q_n|}\)ですが、これは上記の級数を収束させていく事で自然に求まりますね。一応、\(|Q'_n|\)の求め方について以下のように補足します。
\[ |Q'_n| = |Qd_n| = \sqrt{|Qd_n|^2} = \sqrt{ Qd_{n(w)} \cdot Qd_{n(w)} + Qd_{n(x)} \cdot Qd_{n(x)} + Qd_{n(y)} \cdot Qd_{n(y)} + Qd_{n(z)} \cdot Qd_{n(z)}} \]……本当かな？ ここまでわかったら後は以下の式に代入するだけです(再掲)
\[ d \ge \frac{G(\vec{p})}{|G'(\vec{p})|} = \lim_{n \to \infty} \frac{|Q_n| \log|Q_n|}{|Q'_n|} \]\(d\)についてる\(\ge\)が気持ち悪いのですが、実用上は以下のように近似してしまっていいようです。
\[ d \fallingdotseq \lim_{n \to \infty} \frac{1}{2} \frac{|Q_n| \log|Q_n|}{|Q'_n|} \]この1/2という数にまつわる話はソースコードのコメントにこっそり書いてありますので興味があればそちらをどうぞ。
次に、法線の求め方について。今まで描こうとしていた物って結局のところ、\(G(\vec{p})=0\)になる等値面なわけですから、\(\nabla G(\vec{p})\)がそのまま法線になりそうです。さっきは大きさだけ見ていましたが、向きも考慮した\(\nabla G(\vec{p})\)そのものの求め方は
\[ \nabla G(\vec{p}) = ( \frac{\partial G(\vec{p})}{\partial w}, \frac{\partial G(\vec{p})}{\partial x}, \frac{\partial G(\vec{p})}{\partial y}, \frac{\partial G(\vec{p})}{\partial z} ) \]とりあえずx軸方向の偏導関数を求めてみましょうか。合成関数の偏微分は手計算だと面倒ですし、折角なので自動微分的に偏導関数を求めてみましょう。
\[ \lim_{n \to \infty} Q_n = Q = \langle Qr, Qd_{(x)} \rangle と書くことにして、\\ |Q|^2 = Q \cdot Q = \langle Qr \cdot Qr, 2Qr \cdot Qd_{(x)} \rangle \\ |Q| = \sqrt{|Q|^2} = \langle |Qr|, \frac{Qr \cdot Qd_{(x)}}{|Qr|} \rangle \\ G(\vec{p}) = \log|Q| = \langle \log |Qr|, \frac{Qr \cdot Qd_{(x)}}{|Qr|^2} \rangle \\ \therefore \frac{\partial G(\vec{p})}{\partial x} = \frac{Qr \cdot Qd_{(x)}}{|Qr|^2} \]他の軸に関しても軸だけ変えれば全く同じでいけます。法線は最終的に正規化されるので、x,y,z,w方向それぞれを同じスカラー倍にしてる部分は削って良いですから法線を\(\vec{N}\)とすると
\[ \vec{N} = normalize( \nabla G(\vec{p}) ) = normalize(\lim_{n \to \infty} (Qr_n \cdot Qd_{n(w)} , Qr_n \cdot Qd_{n(x)}, Qr_n \cdot Qd_{n(y)}, Qr_n \cdot Qd_{n(z)}) ) \]※GLSLで書くときはvec4は(x,y,z,w)の順になるので注意してください
が、頑張った……法線の求め方に納得が行かなかったのでああでもない、こうでもないとしていたら一週間近く使ってしまいましたが、なんとか求めたかった式にたどり着けました。
わざわざ\(\lim_{n \to \infty}\)が付いてますが、「漸化式の繰り返し計算を終えた後の値ですよ」という事を強調するためこうしているだけの話ですから、別に恐れることは何もないと思います。
今回は画面上の座標(x,y,z)から\(\vec{p} = (z, 定数, x, y)\)として指定しましたから、画面上の出力と矛盾しないように Normal = N.yzw などして3次元の法線に戻す際に軸を入れ替えてやる必要があります。これを忘れるとハマります。大事なことなので二回言いました。
今回の文書の作成には以下を参考にさせていただきました。ありがとうございました！
やっぱりブログを始めてトラックバックとかするようにした方がいいのかなぁと考える昨今です。
基本的にひとりごちはアップロードしたらよほど見過ごせないミスでもない限りは書いたらそのままで放置しておりますが、今回の記事はとても長い上に見づらいのでちょくちょく修正を入れています。
Hayase Taku(SANDMAN)
戻る