2003年08月06日
ポリゴン面の接線ベクトル
異方性BRDFを実装するときに、接線ベクトル(tangent vector)を求める必要が出てきました。
接線ベクトルは、3次元のサーフェスに垂直な一意の平面(接平面)を構成するベクトルです。平面は2つの基底となるベクトルで表現できます。CGの場合、接平面のx軸となるのが接線ベクトル、y軸となるのが従法線ベクトル(binormal vector)と定義されます。接線ベクトルは、3次元のサーフェスを微分したものであるので、物理学ではグラディエント(gradient, 勾配)ベクトルとも呼ばれます。
図の例では、n がサーフェスの法線ベクトル、u が接線ベクトル、 v が従法線ベクトルになります。n, u, v はそれぞれユーグリッド空間で空間を定義する z, x, y ベクトルに対応します(右手座標系)。
これまでも、ライティングの時に半球上にランダムにサンプル点を生成するときに、ローカル座標(z-up)でサンプル点を生成しておいてそれをz軸がサーフェスの法線方向を向くようにサンプル点を座標変換するときに、この接線ベクトルを求める必要がありました。
しかし、このときは、サンプル点を変換するための座標系の基底ベクトルが必要なだけであって、接平面で接線ベクトルがどの方向を向いているかというのは重要ではありませんでした。
しかし、異方性BRDFの場合には、接平面での接線ベクトルの方向をちゃんと求めなければなりません。
サーフェスがパラメトリック曲線であれば、法線ベクトルとともにuvパラメータ値があるので、接線ベクトルと従法線ベクトルは自然と求められますが、lucille ではジオメトリとして三角形のみしか扱わないので、ポリゴン面から接線ベクトルを求める必要があります。
ここで、ポリゴン面の接線ベクトルを求めるときに2つのケースが考えられます。
1. テクスチャ座標がある
2. テクスチャ座標がない
1. のテクスチャ座標があるときは、これをパラメトリック曲線のuvパラメータとして考えることで、接線ベクトルの近似を求めることができます。
2. のテクスチャ座標がない場合は、これはもうメッシュやポイントレンダリングの分野の問題になります。私はあまり詳しくはありませんが、たとえば主成分分析などで求めたりすることができると思います。また、
Anisotropic Polygonal Remeshing
Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Levy and Mathieu Desbrun, SIGGRAP 2003
あたりが参考になるかと思います。
以下では、1. の場合を考えます。
ポリゴンの頂点と、テクスチャ座標から、接線ベクトルを求める方法は、インタラクティブCG分野でポリゴンメッシュからバンプマッピングのための基底ベクトルを求めるのと同じです。
以下が参考になります。
The Differential Geometry of Texture Mapping and Shading
Ken Turkowski, Apple Computer Technical Report, 1992
このレポートでは、ポリゴンメッシュ上のありとあらゆる微分幾何学についての解説がされています(奥がふかいなぁ)。この一部に、ポリゴンメッシュから接線ベクトルを求める方法が掲載されています。
接線ベクトルが求まれば、従法線ベクトルは法線ベクトルと接線ベクトルとの外積で求めることができます。
ポリゴン面から接線ベクトルを求めるのはあくまで近似です。ポリゴンが密であればうまくいくでしょうが、ポリゴンが疎であまにもゆがんだテクスチャ座標が割り当てられていたりすると、当然ですがおかしくなります。
2003年10月24日
数学ソフトウェア
コンピュータグラフィックスは数学のかたまりです。
プログラムを書いているときに、いつも直面するのが数学部分のコードで、
このコードはホントに正しく動作しているのか? 計算が合っているのか?
でいつもものすごい時間を取られます。
このようなコードは、リアルタイムCGであれば行列・ベクトルの計算、物理シミュレーションであれば微分の計算、レンダリングであれば積分の計算が主になるでしょう。
たとえば、マウスの動きからカメラの行列を算出したり、主成分分析で有向バウンディングボックス(OBB)を計算したり、関数のテーラー展開したりなどが挙げられます。
実際にコードが正しい値を出しているかを確認するのは、gdb や printf() でデバッグするのも悪くはありませんが、おなじ計算を数学ソフトウェアで計算させてみるのが一番です。
既存の数学ソフトウェアで検証すれば、偉大なる先人たちのおかげでこれらはまず間違った計算をしないでしょうし、また任意長の精度での計算がサポートされていたりするので、精度的にも信頼できます。
お金があれば微分積分や記号演算は Mathematica、行列・ベクトル計算は Matlab で決まりでしょうが、実際にはそれほどいろいろな機能を多用することもないので、フリーで入手できるソフトで十分な場合がほとんどだと思います。
記号演算など、Mathematica に相当するのは Maxima があります。
http://maxima.sourceforge.net/
行列計算など Matlab に相当するのは gnu octave があります。
octave は簡単な電卓としても使えますし、readline を利用しているので、bash 系を利用しているひとは使いやすいと思います。
ただ、Mathematica は学生であれば学生版が 3 万円くらいで Win, Mac, Linux 版が収録されたハイブリッド版が買えるので、結構お得だとおもいます。
2003年11月22日
クリフォード代数でレイトレ
NICOGRAPH で幅優先レイトレの権威と会うことができました。
で、そのときに教えていただいたもので、おもしろそうなのがこれ。
Modelling 3D Euclidean geometry
by D. Fontijne and L.Dorst, part III
IEEE Computer Graphics and Applications(CGA) March/April 2003.
IEEE 出版物は全然守備範囲外だったなぁ...これからは要チェックですね。
ぱっと見たところ、この論文では、3 次元の幾何を超えて、4 次元、5 次元の世界で幾何を考えていろいろスッキリ(?)記述しようって感じのものです。
幅優先レイトレ(ray-z-bufferベース)では、レイとプリミティブの交差判定などの幾何についての洞察が必要になるのでこのような論文というものはよいアイデアになります。
で、実際にはクリフォード代数(Clifford algebra)や4D プリュッカー(Pluecker)座標でレイトレに必要な幾何を記述するということが書かれています。
クリフォード代数は、geometric algebra としても知られています。geometric algebra の訳語は見当たらないけど幾何(的)代数って感じでしょうか。しかしこれだと代数幾何とかぶってまぎらわしくもありますが。
この論文をはじめとして、クリフォード代数(とそのグラフィックスへの応用)については、以下のサイトにまとめられています。
Geometric algebra(Clifford algebra)
クリフォード代数のライブラリなどが公開されていますのでかなり参考になるかと思います。
論文に記載されているクリフォード代数などを用いたレイトレーサは、以下で入手できます。
http://carol.science.uva.nl/~fontijne/raytracer/
2003年12月01日
作用素とは?
作用素(operator) というのがよく判らなかったのですこし調べてみることにしました。
ってか岩波数学事典を持っているのになんでこれを引くことにに気づかなんだ...
岩波数学事典によると、
(実または複素)線形空間 X の部分集合 D(T) から線形空間 Y への写像 T を X から Y への作用素(operator)と呼ぶ。
とのこと。また MathWorld での作用素の定義。
http://mathworld.wolfram.com/Operator.html
ついでに今回のレンダリング方程式(第二種フレッドホルム方程式)とからむ作用素の歴史を岩波数学事典から
--
関数解析の起源は、1887 年に始められた V. Volterra の研究にさかのぼると考えられる。
彼は、その定義域も地域もともに関数の集合であるような、一般化された関数の概念の重要性を強調した。
このような操作または作用素(operator)の典型的な例は、関数 f にその導関数 f' を対応させる作用素である。
そして、それが退化したものとして、関数 f に数値 f(a) や数値 ∫_{a}^{b} f(t) dt を対応させるものがあり、これは今日 f の汎関数(functional)と呼ばれている。
1896 年に Volterra は、連続関数 f に対して積分方程式 f(x) = ψ(x) - ∫_{a}^{x} K(x,y) ψ(y) dy の連続解 ψ を対応させる作用素を考えた。ここに K(x,y) は x, y の連続関数である。彼は恒等作用素 I (Iψ = I) と (K ψ)(x)= ∫_{a}^{x} K(x,y) ψ(y) dy で与えられる作用素 K を用い、
ψ = (I - K)^{-1} f = f + K f + K^2 f + ... = K^n f = K(K^{n-1} f)
によって解 ψ が求められることを示した。
これに示唆された I. Fredholm は 1900 年に、パラーメータ λ を含む積分方程式
f(x) = ψ(x) - λ ∫_{a}^{x} K(x,y) ψ(y) dy
を取り上げた。
...
--
とのこと。つまりレンダリング方程式(第二種フレッドホルム方程式)を作用素で書いてノイマン展開するってのは積分の分野の人たちには良く知られた有名なことがらだったんですね。作用素のノルムが 1 より小さいときに、これがノイマン展開できるという証明も調べるとありました。
まあ作用素とは基本的には写像のことと捉えていいのかな。また我々が使う微分記号 dy/dx も作用素(微分作用素)になるようです。
ちょっとこれ以上は奥が深くなるのでここでやめておく事にします。
2003年12月03日
レンダリング方程式のノイマン級数展開 間違い
前回書いた レンダリング方程式のノイマン級数展開 ですが、 よくよく Kajiya の論文を最初からよく読んだら、g はパラメータではなくてジオメトリ項(可視関数を含む) でした。
第二種フレッドホルム積分方程式のパラメータの λ のことを指しているのとてっきり勘違いしていました。
また Kajiya の論文では、積分の区間(測度)は半球(立体角)ではなくてサーフェスでした。
2003年12月02日
測度とは?
測度(measure)という考えが積分や確率を考える時に重要になってくるようで、MLT の論文でも measure という単語がよく出てきており、 MLT を理解する上で測度を理解している必要があるようです。
MLT を始めとして、 モンテカルロ法での確率論やLDS などにも測度が出てきたりしています。
このグラフィックスの分野にもちょこちょこ出てくる測度って何だろうということで、今回も岩波数学事典を引いてみました。
--
大まかに云えば、測度はある空間の部分集合の非負関数で完全加法的なものである。ここで完全加法的とは互いに素な集合の列の和の測度が、それぞれの集合の測度の和となることを意味する。この概念は長さ、面積、体積、質量分布、確率分布などの数学的抽象化として導入されたものである。測度論は積分論の基礎であり、この両理論は現代数学、とくに解析学、関数解析学、確率論において重要な役割を演じている。
--
とのこと。
うーん、つまりは測度というのはマイナスでなく、また要素の集まりが重なり合わないもので、長さ、体積、などを一般化して捉えたもの、ということでしょうか。
完全加法的というのは、りんごが 10 個あったとして、その重さは 1 個 1 個別々に測ってこれを 10 個分足しても、 10 個分一度にまとめて測ってもどちらもりんご 10 個の重さは同じで変わらないってことかな。
2003年11月28日
レンダリング方程式のノイマン級数展開
パスレイトレーシング、MLT をはじめとして、多くのレンダリングのアルゴリズムの基礎となっているのが、 Kajiya によるレンダリング方程式です。
James T. Kajiya
The rendering equation
SIGGRAPH 1986, pp 143-150, 1986
MLT の論文を読んでる時に、レンダリング方程式からノイマン級数展開してできる結果の導出がどうもよく理解できていなかったので、復習することにしました。
レンダリング方程式は、こんな形をしています(論文や書籍により細かい表記はまちまちです)。
L(x → w) は、サーフェスの点 x から w 方向に出て行く放射輝度(照明値)、 L(x ← w) はその反対で、点 x に w の方向から入ってくる放射輝度です。基本的には方向は入れ替えても値は同じですが、ここでは明示的に方向を示すために矢印を用いて書いています。
Le(x → w) は点 x から w 方向へと発光する放射輝度です。これには光源や、自己発光する物体などが当てはまります。基本的には光源のサーフェスでのみ非ゼロで、ほかのサーフェスではゼロになります。
積分の区間 Ω は、サーフェスの法線の向き側の半球を示します。
f() は BRDF です。 θ は、サーフェスの法線と w' との角度になります。
dw' は微小立体角です。 これは球座標(θ, φ) で表現すると、
dw' = sinθ dθ dφ
となります。mathematica などの数値計算で球面上の積分を行う場合にはこの球座標での表現を用います。
レンダリング方程式を図示するとこんな感じになります。
つまり、ある点 x から、ある方向 w への照明値というのは、 w 方向へ発光する照明値と、点 x に半球上のすべての方向から入ってくる照明のうち、w 方向へと反射していく光をすべて考慮したものになります。
今回はレンダリング方程式は半球を基本として書いてありますが、これをシーン全体のオブジェクトサーフェスを基本として定式化する場合もあります。MLT の論文ではサーフェスを基本としたほうになっています。
問題
レンダリング方程式は、第 2 種フレッドホルム積分方程式としても捉えることができます。
で、何が問題なのかというと、この積分方程式には、積分したい部分に、積分した結果が含まれていることです。
ここで、 L(x → w) と L(x ← w) は同じなので、同じ L とすると、
L というのは、 L を積分したもの
ということになります。
つまり L を求めるには、 L を積分することになるが、その L 自身も積分の結果であるので、ループしてしまっているのです。
点 A と 点 B の２点での簡単なこの例を説明すると、
ある点 A の照明値を求めるには、点 A から見える点 B の面の照明値が必要になるが、
その点 B の照明値を求めるためにも、点 A の照明値が必要
ということです。
作用素
この式を普通に解こうとすると、左辺の L を右辺の積分内の L へと代入を繰り返していくことになります。
ただ、その場合だと非常に式が長ったらしくなるので、ここで作用素(operator)と呼ばれる考えを導入します。
operator の訳語は、物理学分野では演算子、数学分野では作用素となっているようですが、ここでは数学の方を採用することにします。
今回の作用素 T は、積分の式で L を抜いたものを表します。
・の部分に作用素がかかるものが入ります。 TL とすると、・の部分には L が入ります。まあ関数みたいなものですね。
例えるなら T がヘブンズドアーで、 T のかかる L が書き込む内容という感じでしょうか。
ヘブンズドアー! 「今のことは忘れる」 -> ホントに忘れる
ヘブンズドアー! 「イタリア語をは話せるようになる」 -> ペラペラ
ヘブンズドアー! 「地獄へ行く」 -> 念のため描いておいてやるよ
この作用素の形式でレンダリング方程式を記述すると、
L = Le + TL
という風に簡潔にまとめることができます。ここから、以下では x と w は自明であるとして省略することにします。
ここで、右辺の TL の L に、 L = Le + TL を代入することで、
L = Le + T(Le + TL) = Le + TLe + T^2 L
を導けます。この代入の展開を 2 回行うと、
L = Le + T(Le + T(Le + TL)) = Le + TLe + T^2 Le + T^3 L
になります。これを無限に繰り返していけば、
L = Le + TLe + T^2 Le + ... = Σ_{n=0}^{∞} T^n Le
と書く事ができます。
この式の重要なことはなにかというと、 L を求めるのに、 L に依存せずに Le のみに依存する式へと変形することができるということです。
ノイマン級数展開
で、ここからがよく分っていなかった、レンダリング方程式のノイマン級数(Neumann series)展開です。
結果的には上の L = Σ^{∞} T^n Le と同じになります。
まず、ノイマン級数ってなんなのか知りませんでした。調べたら、ノイマン級数とは、
を云うことでした。いわゆる等比級数の公式だったんですね。ここで 0 <= z <= 1 です。
まず、作用素でのレンダリング方程式の記述に定数項 g を加えて一般的に考えます。
L = gLe + gTL
Kajiya のオリジナルの論文でもそうなっているのでこうしました。
まず、右辺の gTL を左辺に持っていきます。
L - gTL = gLe
左辺を L でまとめると、
(1 - gT)L = gLe
ここで左辺を L だけの式にすると、
L = (1 - gT)^{-1} gLe
になります。ここで作用素も通常の変数と同じに扱えば(ここらへんの数学的な確証はよく分りません)、ノイマン級数の式が適用できるので
が導けます。ノイマン展開では展開する変数が 0 <= z <= 1 でしたが、これが作用素 T にあてはまるかということですが、エネルギー保存則が成り立つのであれば、つまり BRDF積 f() Cosθ の積分が 1 以下になるのであれば、作用素の大きさは 1 以下になるので適用することができると解釈しました。多分これについてはなんらかの証明があるんでしょうね。 まあ証明を知らなくても、1 より大きな定数になった場合は g でスケールして 1 より小さくすればよいし、いずれにせよ積分内ではエネルギー保存則により入力の照明値より出力の照明値が大きくならないと考えれば、なんとなく分る気がします。
また 式を見ても、総和は非負の単調増加なので、 gTLe < gLe < 1 が成り立つひつようがあるので、そこからエネルギー保存則が適用できると考えて導出することもできるかと思います。
この導出により、 L を Le から、つまりライトからのエネルギーのみから、積分計算を行うことで、照明値を計算することができることを示すことができます。
ここから、Kajiya の論文ではパストレーシングの導出とかをしているのですが、それはまた別のお話ということで。
2004年01月27日
放射輝度とは
放射輝度(radiance)を始めとして、放射測定(radiometry)で用いられる単位量やその積分のしかたについて結構納得がいかない部分があり、ずっと考えています。たとえば、放射輝度。放射輝度は、
と定義されています。多くのテキストではこの定義を天下り式に用いているのですが、直感的な疑問が一つ。cosθ の θ は点 x での法線と立体角 ω の間の角になります。これは放射輝度が単位射影面積(cosθdA)で測るため含まれるのですが、cosθは微分作用素から外れているので、ここで θ -> π/2、つまり地平面(法線と直角)の方向とすると、cosθ= 0 となり L が無限大を取ることになってしまいます。
しかし実際には L は有限の値を取ります。なのでどこかで無限大が相殺される部分があるのですが...
L は BRDF積(f cosθ) で積分することでも求められますが、このときは cosθ が打ち消されるのでこのような疑問は生じません。しかし L を上記の式から直接求めたい(光源から発する放射輝度の量など)と思うとこの疑問の壁にぶつかるのです。
同様に立体角の緯度経度表現 dω=sin θ dθ dφ も、これを用いて積分を行う場合、通常のリーマン的な積分(積分区間を短冊に区切って近似する)を行おうとすると、 θ=0 のとき、つまり法線方向の関数の評価を行おうとすると sinθ=0 になるので関数評価がゼロになってしまうのではないかという疑問があります。
たぶん自分の積分微分の考え方が間違っているとおもうので、Stanford の講義ノートや Advanced Global Illumination、 Eric Veach 博士の博士論文などを読んでいるところなのですが、いまいちまだ納得する解に至っていません。積分微分をもう一度やり直さないとダメですね。
2004年02月03日
測度と確率
メトロポリス光輸送(MLT)やパストレーシング(PT)の導出における基礎の理論として、まずレイの経路(path)の測度(measure)という問題にぶつかるかと思います。
今は MLT, PT の実装および基礎の理解として、Eric Veach の博士論文を読んでいます[1]。というかこの博士論文、内容も濃いのに 400 ページ以上もあって凄すぎです。
博士論文によれば、MLT や PT の核となるのは光輸送を経路積分(Path integration)という積分の形で表現するということです。積分の形で表現できれば、既存の積分のテクニックであるモンテカルロ法などを適用することができます。
しかしそのときに考慮しなければいけないことは、レイがたどる経路の積分をどう定義するのか、ということです。
積分を定義するには、積分の区間(domain, 積分記号 ∫ の上や下に付くもの)、積分の測度(measure, dx, dy, df(x) など積分式の最後に付くもので d を取ったものに相当)、被積分関数(integrand) の 3 つが必要になります。
PT, MLT では、積分の区間はレイがたどりうるすべての経路(もしくはシーンのサーフェス全体)になります。被積分関数は BRDF f や 可視関数 V, 放射輝度 L などになります。
そして、ではレイ(の経路)の測度とは? となるのですが、まずは確率における測度という考えから始めた方がわかりやすいかと思います。
測度とは? も簡単な基礎知識として参照
確率と測度
確率という不確かなものを、測度(測度はルベーグが学位論文! で初めて用いて定義した)という決定的な考えで定義しなおしたのがコルモゴルフです。コルモゴルフにより定義された確率論は近代確率論と呼ばれ、確率という不確定的なものを数学的に取り扱ううえで大きく発展しました。
測度的な考えでサイコロの目の出る確率を考えると、以下のようになります。
まずすべての目の出方の組み合わせ(目が取りうるすべての状態)を考えます。サイコロを一回ふる場合、これはすべてで 6 通りあります。
1, 2, 3, 4, 5, 6
サイコロを 2 回ふる場合、これはすべてで 36 通りあります。
11, 12, 13, 14, 15, 16
21, 22, 23, 24, 25, 26
...
61, 62, 63, 64, 65, 66
おのおのの状態は負ではなく、また同時に同じ状態を取ることもありません(完全加法的)。つまり目の出る状態というのは測度として定義することができます。このとき、 目 22 が出るのは、 36 通りのうち 1 通りになります。確率的に云えば 1/36 です。これを測度的に考えれば、 36 (の面積)のうち 1 (の面積)を占めるもの、として解釈することができます(つまり目 22 は 1/36 の面積)。
実際にはすべての状態を 1 にスケーリングしたもので考えます(確率 1 と定義できるため)。
このようにして、レイのたどるすべての組み合わせ(状態)から、あるレイ(1 本でも数本でも)のすべての状態に占める割合(面積つまりは測度)を考えることができ、それが定義できれば積分を定義することができます。
これについてはまた次回。
[1] Eric Veach, Robust Monte Carlo Methods for Light Transport Simulation Ph.D. dissertation, Stanford University, December 1997.
2004年02月18日
測度論的放射測定
Eric Veach の論文では、いたるところに測度(measure)という考え方による説明が見受けられます。
Veach の博士論文[1]を読んでいて、第三章の付録に測度論的放射測定(Measure-theoretic radiometry)についての説明がありました。以下は三章付録の最初の部分の日本語訳。
通常、放射測定量(radiometric quantities)は"無限小(infinitesimals)"を用いて極限(limit)を取ることで定義されています。 Arvo[2] の 2 章では光子(photon)の観察可能(observable)なふるまいに関連した公理(axioms)を提案することで、これとは異なったアプローチを取っています。つまり測度論の考えを用いて放射測定量を導出しています。 彼の解析は定常状態(steady-state)でモノクロの放射(radiation)での空間分布に焦点を当てています。そして位相空間密度(phase space density)の測度論的な定義を導きだしています。本章では、彼のテクニックをより一般的な放射測定量の種類にも拡張します。例えばスペクトル放射輝度(spectral radiance)やスペクトル放射ステリゼント(spectral radiant sterisent)(訳注: sterisent という単語は何か分りませんでしたが、radiant sterisent については Nicodemus の文献に詳しい)などへの測度論的な定義を与えます。
放射輝度などの定義についての、測度論によるアプローチは James Arvo の博士論文 [2] が元ネタのようです。
多くの CG のテキストでは、放射測定を説明するときには、微分的な考えか、光子の密度を用いて説明しているかと思います。
前回の 測度と確率 で取り上げたように、密度などの確率的なものは、測度を用いてより数学的に厳密に説明することができます。
その意味では、放射輝度なども密度であるので、放射測定を測度論を用いて説明するというのは自然ですね。
いわゆる微分的な考え(無限小を考える)での放射測定の説明だと、僕のように、放射輝度の式の分母に cosθがあると水平になったときに無限大を取るのではないか? とか疑問がでてきてしまったり、ではこの微分の式というのはそのまま評価してもよいのか(放射輝度などは無限小の微分の形そのままでは評価できず、実際には有限の平均でしか観測ないことは Nicodemus が指摘しているが、多くのテキストではそのことには触れていない)などといろいろ不自然にしっくりこないことが多いので、測度論的な放射測定の説明のほうがもし分りやすく説明できのであれば、そちらのほうが変な疑問が湧いたりせずにロバストで落とし穴なく解釈できるのではないかと思います。
Eric Veach が測度を使って解説するのを愛してやまない? のもそこらへんに理由があるのかもしれません。
レンダリングアルゴリズム(とくに光輸送)を数学的に厳密に記述する時に、測度論という考えが非常に重要になるかと思います。測度論は奥が深いので、きっと極めるといい事があるかもしれません。
ところで James Arvo って結構古参な人だと思っていたのですが、博士論文は 95 年と最近なのですね。
[1] Eric Veach, "Robust Monte Carlo Methods for Light Transport Simulation", Ph.D. dissertation, Stanford University, December 1997.
[2] James Arvo, "Analytic Methods for Simulated Light Transport", Ph.D. Thesis, Yale University, December, 1995.
2004年03月21日
rover と線形合同法
MLT では、画像面のピクセルをランダムにサンプルするプロセス(レンズサブパス変異)があります。
これを、ただ単に擬似乱数でサンプルしようとしても、
すべての画像面のピクセルを一様にまんべんなく取り上げていくことはなく、
サンプルされない領域とサンプルされすぎる領域に分かれてしまい、
非常に効率が悪いというのは良く知られています。
この問題には、層別化や LDS を使うというのも考えられていますが、
Eric Veach 氏の博士論文では、これに rover(ローバー)と呼ばれる
メモリ管理法からの応用を利用することを述べられています。
(11.4.4 節)
rover というのは、各ピクセルが一度だけ表れるようなランダムなピクセル順を
生成するインデックスアルゴリズムのことのようです。
(rover memory とかで検索すると火星探査機の rover がヒットしてしまう)
各ピクセルにはサンプルされた回数を記録するカウンタを用意しておき、
カウンタがある値を超えたら(サンプルされすぎたら)、
rover によるインデックスを参照して次のランダムのピクセルを
採用するようにします。
rover では、各ピクセルが一度だけ表れるようなランダムなピクセル順を生成するのに、
線形合同法(linear congruential) の低位ビットを利用しています。
各ピクセルが一度だけ表れるというのは linear pixel shuffling でも同じですね。
(linear pixel shuffling と線形合同法は式も似ていますしね)
線形合同法
線形合同法は、擬似乱数を生成するアルゴリズムの古典的な方法で、たとえば
gcc では以下のように実装されているらしいです。
(たぶん最近では Mersenne Twister に代わっているでしょうか?)
unsigned int seed = 1;
int
randLC()
{
seed = 1103515245 * seed + 12345;
return seed;
}
という非常に単純なアルゴリズムです。
線形合同法では、低位 n ビットを取り出すと、それが周期 2^n の
数列を成すことが知られています。
またこのとき、各周期 2^n では、おのおの重なる数字は現れません。
この性質を利用する事で、各要素の数字が必ず 1 度しか現れないランダムな
数列を生成することが可能になります。
ためしに 2^8 = 256 個の数列を作って重なりがないかテストするコードを書いてみました。
/* lc.c */
#include <stdio.h>
#define GENSIZE 256
unsigned int seed = 1;
int table[GENSIZE];
int
randLC()
{
seed = 1103515245 * seed + 12345;
return seed;
}
void
init()
{
int i;
for (i = 0; i < GENSIZE; i++) {
table[i] = 0;
}
}
void
test()
{
int i;
int r;
int t;
int count = 0;
for (i = 0; i < GENSIZE; i++) {
r = randLC();
t = r & (GENSIZE-1); /* extract lower bits */
if (table[t] != 0) {
printf("boo! t = %d\n", t);
}
table[t] = 1;
}
for (i = 0; i < GENSIZE; i++) {
if (table[i] == 1) count++;
}
printf("count = %d\n", count);
}
int
main()
{
init();
test();
}
実行結果。
$ gcc lc.c $ ./a.out count = 256
と、重なりが無い事が確認できます。
Eric Veach の博士論文では、さらに各ピクセルの内部(サブピクセル)では、
Poisson minimum-disk パターンというのを用いているそうです。
2004年06月15日
non-Negative Matrix Factorization
今年の SIGGRAPH 2004 の論文、Efficient BRDF Importance Sampling Using A Factored Representation では、データを圧縮、インポータンスサンプルしやすいように、二次元の行列を 2 つのより要素数の少ない二次元の行列の積に分解( Y = GF のように、 たとえば 4x4 行列を、 4x1 と 1x4 の 2 つの行列に分解 ) しているのですが、このとき非負値行列分解(non-negative matrix factorization, NMF)と呼ばれる手法を用いています。
NMF は、名前の通り負値の要素のない行列を、分解後の行列もまた負値の要素を持たないようする手法です。特異値分解(singular value decomposition, SVD)や主成分分析(principal component analysis, PCA)なども同じように行列を分解しますが、NMF は負値を含まないので、確率やインポータンスのように負値を取らない対象を扱うときに有益です。
NMF のアルゴリズムは Lee [1] で詳しく述べられています。最近提案された新しいアルゴリズムのようです。NMF は、SIGGRAPH 2002 の論文であるライトフィールドマッピング [2] でも用いられています。
アルゴリズム
詳細な NMF 分解アルゴリズムは [1] に述べられていますが、ここでは簡単に述べたいと思います。
NMF は、反復的なアルゴリズムです。
たとえば Y = GF という風に行列 Y を G と F に分解するとします。
まず、G と F の要素は、非負のランダム値で埋めておきます。そして、あとは更新のルールに従って、 G と F を反復的にアップデートするというものです。処理自体は単純そう。
シェーダとの組み合わせ
Efficient BRDF Importance Sampling Using A Factored Representation の手法は、RenderMan シェーダのインポータンスサンプリングにも使えそうです。しかし、その場合は解決しなければならない問題として、
o シェーダコードを言語解析して、BRDF に相当する部分を抜き出す。
o シェーダのサンプリングや、NMF 分解の時間がどれくらいかかるかを見極めておく必要がある(ただし一回計算できれば OK なので、結果をキャッシュして保存しておくことができる)
が挙げられます。
PRMan 11 では、フォトンマッピング時に、サーフェスの反射のプロパティをいくつかのプリセットから指定できるようです。
Attribute "photon" "string shadingmodel" ["matte"] (glass,water,chrome,matte,transparent,"")
もしプリセットが指定されていない場合は、シェーダコードを解析して反射のプロパティを抽出するようですが、どのようなアルゴリズムが使われているのかは判りません。
グローバルイルミネーションでは、インポータンスが必須です。
既存のシェーダ言語にはこれらインポータンスを記述する機構はありません。シェーディングのコードブロックからシェーダコンパイラがインポータンスを解析して抽出するのは難しいかもしれないので、これからのグローバルイルミネーションを視野に入れたシェーダ言語設計では、ユーザが明示的にインポータンスのコードブロックも記述できるようであったほうがいいのかもしれません。
[1] LEE, D. D., AND SEUNG, H. S. 2000. Algorithms for non-negative matrix factorization. In NIPS, 556–562.
[2] CHEN, W.-C., BOUGUET, J.-Y., CHU, M. H., AND GRZESZCZUK,R. 2002. Light field mapping: efficient representation and hardware rendering of surface light fields. In SIGGRAPH 02, 447–456.
NMF コード
というわけで、NMF を実装してみました。assert() で非負値とゼロ除算のチェックをしていますが、行列のサイズが大きくなるとほとんどゼロになる要素も出てくるので、ここらへんは改善したほうがよいでしょう。 NMF の分解はとても早く、数百 X 数百の行列のサイズを 10 回反復で分解させても、1 秒程度で計算ができます。
/*
* Algorithm implementation of Non-Negative Matrix Factorization.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "mt19937ar-cok.c" /* Mersenne Twister random number generator */
#define randomMT genrand_real2
void
mat_mul(double **out, double **m1, double **m2, int n, int m, int r)
{
int i, j, k;
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
out[j][i] = 0.0;
for (k = 0; k < r; k++) {
out[j][i] += m1[k][i] * m2[j][k];
}
}
}
}
/* measures NMF error. */
void
nmf_error(double **Y, double **GF, int n, int m)
{
int i, j;
double err = 0.0;
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
assert(GF[j][i] >= 1.0e-12);
err += Y[j][i] * log(Y[j][i]/GF[j][i])
- Y[j][i]
+ GF[j][i];
}
}
printf("NMF err = %f\n", err);
}
/*
* Given a non-negative matrix Y(n x m),
* find non-negative matrix factors G(n x r) and F(r x m) such that:
* Y = G F.
*/
void
nmf(double **Y, double **G, double **F, int n, int m, int r)
{
int i, j, k;
int loop;
const int maxloop = 10;
double **GF; /* approximation of the target matrix */
double sum;
double sumGj;
GF = (double **)malloc(sizeof(double *) * m);
for (i = 0; i < m; i++) {
GF[i] = (double *)malloc(sizeof(double) * n);
}
/* initialize G and F with random seeds */
for (j = 0; j < r; j++) {
for (i = 0; i < n; i++) {
G[j][i] = randomMT();
}
}
for (j = 0; j < m; j++) {
for (i = 0; i < r; i++) {
F[j][i] = randomMT();
}
}
mat_mul(GF, G, F, n, m, r);
for (loop = 0; loop < maxloop; loop++) {
printf("iteration %d of %d\n", loop, maxloop);
printf("updating Fij...\n");
for (j = 0; j < m; j++) {
for (i = 0; i < r; i++) {
sum = 0.0;
for (k = 0; k < n; k++) {
assert(GF[j][j] >= 1.0e-12);
sum += G[i][k] * (Y[j][k] / GF[j][k]);
}
F[j][i] = F[j][i] * sum;
}
}
printf("updating Gij 1...\n");
for (j = 0; j < r; j++) {
for (i = 0; i < n; i++) {
sum = 0.0;
for (k = 0; k < m; k++) {
assert(GF[k][i] >= 1.0e-12);
sum += (Y[k][i] / GF[k][i]) * F[k][j];
}
G[j][i] = G[j][i] * sum;
}
}
printf("updating Gij 2...\n");
for (j = 0; j < r; j++) {
sumGj = 0.0;
for (k = 0; k < n; k++) {
sumGj += G[j][k];
}
assert(sumGj >= 1.0e-12);
for (i = 0; i < n; i++) {
G[j][i] = G[j][i] / sumGj;
}
}
mat_mul(GF, G, F, n, m, r);
nmf_error(Y, GF, n, m);
}
for (i = 0; i < n; i++) {
free(GF[i]);
}
free(GF);
}
int
main(int argc, char **argv)
{
int i, j;
int n = 4, m = 4, r = 2;
double **Y, **G, **F, **GF;
Y = (double **)malloc(sizeof(double *) * m);
G = (double **)malloc(sizeof(double *) * r);
F = (double **)malloc(sizeof(double *) * m);
GF = (double **)malloc(sizeof(double *) * m);
for (j = 0; j < m; j++) {
Y[j] = (double *)malloc(sizeof(double) * n);
F[j] = (double *)malloc(sizeof(double) * r);
GF[j] = (double *)malloc(sizeof(double) * n);
}
for (j = 0; j < r; j++) {
G[j] = (double *)malloc(sizeof(double) * n);
}
for (j = 0; j < r; j++) {
for (i = 0; i < n; i++) {
G[j][i] = 0;
}
}
for (j = 0; j < m; j++) {
for (i = 0; i < r; i++) {
F[j][i] = 0;
}
}
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
Y[j][i] = j * 8 + i + 1;
}
}
nmf(Y, G, F, n, m, r);
mat_mul(GF, G, F, n, m, r);
printf("\nG = \n---\n");
for (j = 0; j < r; j++) {
for (i = 0; i < n; i++) {
printf("%f ", G[j][i]);
}
printf("\n");
}
printf("---\n");
printf("\nF = \n---\n");
for (j = 0; j < m; j++) {
for (i = 0; i < r; i++) {
printf("%f ", F[j][i]);
}
printf("\n");
}
printf("---\n");
printf("\nY = \n---\n");
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
printf("%f ", Y[j][i]);
}
printf("\n");
}
printf("---\n");
printf("\nGF = \n---\n");
for (j = 0; j < m; j++) {
for (i = 0; i < n; i++) {
printf("%f ", GF[j][i]);
}
printf("\n");
}
printf("---\n");
}
2004年06月18日
Safe and effective importance sampling
QMC や RQMC(Randomized QMC) で有名な Art B. Owen による、インポータンスサンプリング(importance sampling, 加重サンプリング、重点的サンプリング、また統計分野では要点抽出法などとも訳される)の解説を統括的にまとめてある論文。ただ単にインポータンスサンプリングを適用しただけでは、逆に分散を増やしてしまうことがあるため、それらを解決する近年発案されたグラフィックスでも有名な Veach の多重インポータンスサンプリング(multiple importance sampling)などをサーベイし、インポータンスサンプリングへの深い洞察を与えています。また、多重インポータンスサンプリングに修正を加えることで、負値要素を持つ積分にもインポータンスサンプリングが適用できることを示しています。インポータンスサンプリング一つ取っても、いろいろ考えることはありそうだと思いました。
Art B. Owen and Yu Zhou
Safe and effective importance sampling
Technical Report, Stanford University, 1999
論文 Abstract 日本語訳
インポータンスサンプリングは、シミュレーションにおいて分散(variance)を減らすために広く用いられている手法です。
しかしながら、安直なインポータンスサンプリングでは逆に大きく分散を増やしてしまうことがあります。
本論文は、この問題を解決するために近年発案された手法である、Hesterberg の defensive mixture sampling (1995) と Veach の多重インポータンスサンプリング(multiple importance sampling) (1997) をサーベイし、拡張します。
我々が提案する手法は、複合要素の比率を複合密度へと制御変量(control variates)を行ない、複合密度からインポータンスサンプリングを行なうことです。
提案する手法は、複合要素の最適なものからのインポータンスサンプリングよりも悪くなることはなく、より良くなる場合があります。
インポータンスサンプリングは、非負値積分に対してはほぼ完全に適用できることがよく知られています。我々は、多重インポータンスサンプリングに修正を加えることで、これが正と負の値の両方を取る積分にも当てはまることを示します。
2004年07月03日
デカルト座標系での球面調和関数
放射輝度キャッシュ(radiance caching)でまた、もてはやされてきた球面調和関数(spherical harmonics)。
これらの計算には通常は GSL(Gnu Scientific Library) を用いて(θ, Φ)座標で計算しても良いのですが、低次元までであれば球面調和関数は(x, y, z)のベクトルから直接多項式表現で計算した方が速いのでインタラクティブグラフィックス系では後者が好まれるでしょう。
放射輝度キャッシュでも使われるのは l=3,4 ぐらいまでで十分のようですし。
というわけで 1 年くらい前のとおんなじですが、実装に参考になる球面調和関数の多項式表現のリスト(l<=4) までです。
途中までの式は間違ってる可能性もありますが、デカルト座標表現の式自体はたぶん間違っていないはず...