Tweet
Tweet

Slide 1

Slide 1 text

双方向パストレーシングレンダラ edubpt解説 2013/11/19 Ver.1.00 hole(@h013) http://kagamin.net/hole/edubpt/index.htm

Slide 2

Slide 2 text

No content

Slide 3

Slide 3 text

edubptとは シンプルでコンパクトな 双方向パストレーシングレンダラ

Slide 4

Slide 4 text

edubptとは 物理ベース・パストレーシングレンダラ eduptの続編 http://kagamin.net/hole/edupt/index.htm https://github.com/githole/edupt

Slide 5

Slide 5 text

edubptとは C++で記述 教育的 日本語によるコメント付き eduptをベースに双方向化 主な特徴

Slide 6

Slide 6 text

edubptとは Githubでコード公開 https://github.com/githole/edubpt

Slide 7

Slide 7 text

edubptとは 細かい仕様 双方向パストレーシングによる物理ベースレンダラ 物理カメラモデルの採用による被写界深度の再現 リファレンス用のパストレーシングも内蔵（eduptとほぼ同じ） パストレーシング・ライトトレーシングも内蔵（双方向パストレ用） OpenMPによる並列化 標準ライブラリ以外の依存ライブラリは（基本的に）無し

Slide 8

Slide 8 text

このスライドの内容 edubptのコード解説 ＋ 双方向パストレーシングの解説 ＋ その他、物理ベースレンダリングの発展的解説 ※前回の「物理ベースレンダラedupt解説」の内容をある程度理解している前提の内容 ※割と数学の話が多め

Slide 9

Slide 9 text

目次 1. レンダリングに関する基礎知識 1. ロシアンルーレット 2. レンダリングのための確率論 3. 確率密度の変換 4. ジオメトリファクタ 5. マルチプルインポータンスサンプリング 2. レンダリングの定式化 3. パストレーシング 4. 双方向パストレーシング・理論編 1. 双方向パストレが解く式 2. レンダリング方程式 5. 双方向パストレーシング・実装編その１ 1. パストレーシングによるパスのサンプリング 2. ライトトレーシングによるパスのサンプリング 6. 双方向パストレーシング・実装編その２ 7. レンダリング結果

Slide 10

Slide 10 text

1. レンダリングに関する基礎知識

Slide 11

Slide 11 text

レンダリングに関する基礎知識  レンダリングという問題を解いていく際に、必要となる各種の基礎的な理 論・技法について説明していく。 – この章で説明した事柄は、このスライド全体にわたって使う。

Slide 12

Slide 12 text

1.1 ロシアンルーレット

Slide 13

Slide 13 text

ロシアンルーレット  以下のような無限級数を計算したいとする。 – （レンダリングは結局無限級数の計算に帰着される）  計算機は無限級数を扱うことが出来ない。 – 無限ループ  そこで、確率的に処理を打ち切ることで統計的な期待値は真値になるよう にしつつ、現実的な計算量で推定値を得ることを考える。 𝐶 = 𝐶0 + 𝐶1 + 𝐶2 + ⋯ 𝐶 ≈ 𝐶 = 推定値 ただし、𝐸 𝐶 = 𝐶

Slide 14

Slide 14 text

ロシアンルーレット  次のような項を考えると、  𝐶を以下のように書き換えることが出来る。 𝐷𝑖 = 𝐶𝑖 + 𝐶𝑖+1 + 𝐶𝑖+2 + ⋯ 𝐶 = 𝐷0 = 𝐶0 + 𝐷1 = 𝐶0 + 𝐶1 + 𝐷2 = ⋯

Slide 15

Slide 15 text

ロシアンルーレット  以下のように処理を分岐させて𝐶を推定する。 – 𝑟0 は [0, 1)の範囲の乱数 – 𝑝0 は [0, 1]の範囲の定数 – 𝐶 は𝐶の推定器で、確率変数  すると、 𝐸 𝐶 = 𝐷0 = 𝐶となる – 実際の実装では、乱数𝑟0 を何らかの方法で生成し、 𝑝0 未満なら𝐷0 の評価を行い、 𝑝0 以上なら処理を打ち切って𝐶 の評価値を0とする。 – 𝑝0 は基本的に任意の値で良いが、 𝐷0 の値の大小に比例させたほうが分散𝑉[𝐶 ]は小 さくなる。 𝐶 = 𝐷0 𝑝0 , 𝑟0 < 𝑝0 0, 𝑟0 ≥ 𝑝0

Slide 16

Slide 16 text

ロシアンルーレット  𝐷0 の評価を行うことになったとする。  𝐷0 も𝐶 と同じく無限級数なので、再び確率的に処理を分岐させて推定する。  さっきの式と合わせると、推定のための式全体は以下のようになる。 – それぞれの条件分岐のたびに乱数生成を行い、処理を打ち切るか、次の項を評価するかを決める。 – 以下、𝑫𝒊 の評価を行うたびに、それ以降の項𝑫𝒊+𝟏 を評価するかどうかを確率的に決める。 – 𝐸 𝐶 = 𝐸 𝐷0 = 𝐶0 + 𝐸 𝐷1 = 𝐶となる。（以降も同じで、常に推定式全体の期待値は𝐶 になる） – さっきと同様、𝑝𝑖 は基本的に任意の値で良いが、 𝐷𝑖 の値の大小に比例させたほうが分散は小さくなる。 𝐷0 = 𝐶0 + 𝐷1 𝑝1 , 𝑟1 < 𝑝1 0, 𝑟1 ≥ 𝑝1 𝐶 = 𝐷0 𝑝0 , 𝑟0 < 𝑝0 0, 𝑟0 ≥ 𝑝0 = 𝐶0 + 𝐷1 𝑝1 , 𝑟1 < 𝑝1 0, 𝑟1 ≥ 𝑝1 𝑝0 , 𝑟0 < 𝑝0 0, 𝑟0 ≥ 𝑝0

Slide 17

Slide 17 text

ロシアンルーレット  まとめると、推定は以下のように行われる。 𝐶 ≈ 𝐶 = 𝐷0 𝑝0 , 𝑟0 < 𝑝0 0, 𝑟0 ≥ 𝑝0 𝐷𝑖 = 𝐶𝑖 + 𝐷𝑖+1 𝑝𝑖+1 , 𝑟𝑖+1 < 𝑝𝑖+1 0, 𝑟𝑖+1 ≥ 𝑝𝑖+1

Slide 18

Slide 18 text

ロシアンルーレット  𝐶を推定するための疑似コード 1. 𝐶 ← 0, 𝑖 ← 0, 𝑝𝑡 ← 1 2. 𝑤ℎ𝑖𝑙𝑒 𝑡𝑟𝑢𝑒 3. 𝑝𝑡 ← 𝑝𝑡 ⋅ 𝑝𝑖 4. 𝑟𝑖 ← 𝑟𝑎𝑛𝑑𝑜𝑚([0,1)) 5. 𝑖𝑓 𝑟𝑖 < 𝑝𝑖 𝑡ℎ𝑒𝑛 6. 𝐶 ← 𝐶 + 𝐶𝑖 𝑝𝑡 7. 𝑖 ← 𝑖 + 1 8. 𝑒𝑙𝑠𝑒 𝑏𝑟𝑒𝑎𝑘 9. 𝑒𝑛𝑑 10.𝑟𝑒𝑡𝑢𝑟𝑛 𝐶

Slide 19

Slide 19 text

ロシアンルーレット  以上のように、確率的に処理を打ち切るかどうかを選択することで、トー タルの計算量を現実的な量に抑えつつ、推定値𝑪 の期待値は真値と等しくす る（Unbiased）ことが出来る。 – 例えば、𝑝𝑖 = 1 2 とすると、100回連続で処理が行われる（打ち切られない）確率 は 1 2 100 となる。確率的には、100兆回、100京回といった、現実的でない回数連 続で処理が行われてしまう確率も0ではないが、極端に低い確率になるため、そ のようなことは起こらないとしてよい。  レンダリングの文脈においては、シーンをレイトレーシングしてパス（光 路）を追跡していく際、反射回数が無限に増えてしまうのを防ぐためにロ シアンルーレットによって処理を打ち切る。 – また、フレネル反射において反射と屈折のどちらを追跡するかを決定するのも ロシアンルーレット。

Slide 20

Slide 20 text

ロシアンルーレット  今回の例では全ての項の計算についてロシアンルーレットを行ったが、最 初の数項については𝑝𝑖 を1に設定することで打ち切らないようにする、と いったことを行うと最初の数項については、ロシアンルーレットによる分 散が無くなり、トータルの分散を小さく抑えることも出来る。eduptでは最 初の数回のトレースではロシアンルーレットを行っていない。  ただし、edubptにおいては頂点ごとに完全に独立に𝑝𝑖 を設定したいので、 このようなことは行っていない。

Slide 21

Slide 21 text

1.2 レンダリングのための確率論

Slide 22

Slide 22 text

レンダリングのための確率論  確率分布（あるいは確率測度）について考える。  確率分布𝑃は以下のような測度関数として定義される。  𝑃 𝐷 = 𝑃𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦{𝑋 ∈ D}（ここで𝑃はΩ上に定義される測度で、右辺は確 率変数𝑋が集合𝐷に含まれる確率）  対応する確率密度関数𝑝はΩ上の測度𝜇を使ってRadon-Nikodym微分を行うと 𝑝 𝑥 = 𝑑𝑃 𝑑𝜇 (𝑥)となり、これは𝑃 𝐷 = 𝑝 𝑥 𝑑𝜇 𝑥 𝐷 (= 𝑑𝑃 𝑑𝜇 (𝑥)𝑑𝜇 𝑥 𝐷 )を満た すものである。  一般に、測度𝝁に関する確率密度関数𝒑 𝒙 を𝑷𝝁 (𝒙)と表記することにする。

Slide 23

Slide 23 text

具体例  シーンのジオメトリ表面全体を𝑀とする。  𝑀上の確率分布𝑃を考える。  𝑀上の面積測度𝐴に関する確率密度関数を𝑃𝐴 (𝑥)と書く。 – 面積測度とは、ある領域を与えると（いわゆる普通の）面積を返す関数 – 𝑃𝐴 (𝑥)は𝑀上のある一点𝑥についての面積測度に関する確率密度関数となる。 – 今後もしばしば出てくる。  確率分布と測度に応じて様々な確率密度関数が考えられる。 – 同じ面積測度𝐴を使っても、別の確率分布を与えれば、別の確率密度関数になる。 – 𝑃𝐴 (𝑥)と𝑄𝐴 (𝑥)といったように、別の確率密度関数が得られる。 – 同じ確率分布を使っても、別の測度を与えれば、別の確率密度関数になる。 𝑀 𝐷 𝐴 𝐷 = 𝐷の面積

Slide 24

Slide 24 text

具体例2  単位球面全体を𝑆2とする。  𝑆2上の確率分布𝑃を考える。  立体角測度𝜎に関する確率密度関数を𝑃𝜎 (𝜔)と書く。 – 立体角測度とは、単位球面上のある領域を与えると、その領域の立体角を返す関数 – 𝜎(𝐷)は立体角なので単位はステラジアン。 – 𝑷𝝈 (𝝎)は𝑺𝟐上のある一点についての確率密度を表すが、単位球面上の一点というのはすなわ ちある一方向のことなので、 𝑷𝝈 (𝝎)はある一方向についての確率密度関数になる。  確率密度関数を考えるとき、使う測度は（ 𝑆2上の適切な測度なら）なんで もいい。 – 普通の立体角測度𝜎や、投影立体角測度𝜎⊥を使うことができる。 – このスライドでは主に立体角測度を使うことにする。 𝜎(𝐷) = 𝐷の立体角（単位球面上の面積） 𝐷 𝜎⊥(𝐷) = 𝐷の投影立体角（立体角を投影した面積）

Slide 25

Slide 25 text

具体例2  全球から一方向を一様にサンプリングするとき（よくあるシチュエーション）  一様にサンプリングということは確率密度が全方向で定数なので  確率分布の定義より𝑃 𝑆2 = 1となるため  より、 𝑪 = 𝟏 𝟒𝝅 となり、 𝑷𝝈 𝝎 = 𝟏 𝟒𝝅 となる（単位球面の面積の逆数になる）。 – ここで、微小立体角𝑑𝜎 𝜔 が球面上の微小領域sin𝜃𝑑𝜃𝑑𝜙と等しいことを利用して変数 変換して上の式を解いた。 𝑃𝜎 𝜔 = 𝐶 𝑃𝜎 𝜔 𝑆2 𝑑𝜎 𝜔 = 𝐶 𝑆2 𝑑𝜎 𝜔 = 𝐶 𝑑𝜎 𝜔 𝑆2 = 4𝜋𝐶 = 1 𝑑𝜎 𝜔 𝑆2 = sin 𝜃 𝑑𝜃𝑑𝜙 𝜋 0 2𝜋 0

Slide 26

Slide 26 text

具体例2  全球から一方向を一様にサンプリングするとき（よくあるシチュエーション）  同じ分布について、別の測度で確率密度関数を定義してみる。  投影立体角測度𝜎⊥を使う。 – 微小投影立体角𝑑𝜎⊥と微小立体角𝑑𝜎の関係式、 cos𝜃𝑑𝜎 = 𝑑𝜎⊥より  が得られる。𝑃𝜎 𝜔 = 1 4𝜋 であったため上式に代入すると、  となる。同じ一様分布だが、測度によって異なる確率密度が得られるというこ とがわかる。 𝑃𝜎 𝜔 = cos 𝜃 𝑃𝜎⊥ 𝜔 𝑃𝜎⊥ 𝜔 = 1 4𝜋 cos 𝜃

Slide 27

Slide 27 text

モンテカルロ積分  Ωを積分範囲、 𝜇をΩ上に定義される測度とすると、 𝜇に関する確率密度関 数𝑃 𝜇 を使って以下のモンテカルロ積分𝐼 によって上の積分を近似出来る – 𝑋𝑖 は関数に対するサンプルで、 𝑃 𝜇 に従ってサンプリングされる。 𝐼 ≈ 𝐼 = 1 𝑁 𝑓(𝑋𝑖 ) 𝑃 𝜇 (𝑋𝑖 ) 𝑁 𝑖=1 𝐼 = 𝑓 𝑥 𝑑𝜇(𝑥) Ω

Slide 28

Slide 28 text

証明 𝐸 𝐼 = 𝐸 1 𝑁 𝑓 𝑋𝑖 𝑃 𝜇 𝑋𝑖 𝑁 𝑖=1 = 1 𝑁 𝑓 𝑥 𝑃 𝜇 𝑥 𝑃 𝜇 𝑥 𝑑𝜇(𝑥) Ω 𝑁 𝑖=1 = 𝑓 𝑥 𝑃 𝜇 𝑥 𝑃 𝜇 𝑥 𝑑𝜇(𝑥) Ω = 𝑓(𝑥)𝑑𝜇(𝑥) Ω = 𝐼

Slide 29

Slide 29 text

注意  𝐼 = 𝑓 𝑥 𝑑𝜇(𝑥) Ω という積分をモンテカルロ積分で解くた めにはサンプルを測度𝜇に関する確率密度関数𝑃 𝜇 に基づい てサンプリングする必要がある – 別の測度𝜇′に関する確率密度関数𝑃 𝜇′ に基づいてサンプリングした場合、得 られる結果はあくまで 𝑓 𝑥 𝑑𝜇′(𝑥) Ω の近似になり、 𝑓 𝑥 𝑑𝜇(𝑥) Ω の近似に ならない。

Slide 30

Slide 30 text

具体例  単位半球面Ω上に定義されたある関数𝑓 𝜔 を半球面上で微小立体角に関し て積分したいとする。  𝑓 𝜔 𝑑𝜎(𝜔) Ω を解く。  モンテカルロ積分を使うなら、何らかの確率密度関数𝑃𝜎 に基づいてサンプ ルを生成する。  この𝑃𝜎 は、 𝑃𝜎 𝜔 𝑑𝜎(𝜔) Ω = 1を満たしさえすれば、任意の関数を与えるこ とができる（ 𝑓 𝜔 の形状に近い方が誤差の収束が速い→Importance Sampling）  たとえば、 cos 𝜃 𝜋 や 1 2𝜋 などが、条件を満たす確率密度関数となる。

Slide 31

Slide 31 text

1.3 確率密度の変換

Slide 32

Slide 32 text

確率密度の変換 𝑥 𝜔 𝑦  位置𝒙から次の頂点をサンプリングするとする。  レイを飛ばす方向をサンプリングし、その方向にレイトレーシングして、シー ンとの交点が次の頂点𝑦  方向のサンプリングなので、立体角測度に関する確率密度関数に基づいてサン プリングを行う。  例えば、𝑃𝜔 𝜔 = cos 𝜃 𝜋 のような確率密度関数に基づいてサンプリングすれば、cos 𝜃 項を考慮したインポータンスサンプリングということになる。（いつものやつ）

Slide 33

Slide 33 text

確率密度の変換 𝑥 𝜔 𝑦  今回、最終的に解く積分の積分範囲がシーン表面になり、変数（＝頂点）はシーン表面 を動くことになる。（後の章で説明。双方向パストレーシングをすっきり解くためにそ ういう式に変形する）  そのような積分をモンテカルロ積分する際、各頂点は面積測度に関する確率密度関数で サンプリングされなければならないが、実際の実装では立体角測度に関する確率密度関 数によってサンプリングされる。  頂点を直接、面積測度に関する確率密度関数でサンプリングした場合、頂点間が別の物体でさ えぎられて有効なサンプルにならない確率が非常に高くなるため、実際の実装ではレイトレに 基づく立体角測度によるサンプリングを行ったほうが効率が良いため。  よって、立体角測度に関する確率密度を面積測度に関する確率密度に変換する必要が出 てくる。

Slide 34

Slide 34 text

確率密度の変換  位置𝑥における𝜔方向の微小立体角と、対応するシーン内の微小面積との関係式を 立てる。  立体角に対応する面積とは、立体角をシーンに投影したときの、シーン内の領域  図のような幾何的関係が成り立つので、以下の式が得られる。ここで、立体角は単位 球面上の面積に等しいことを利用した。 𝑥 𝑑𝐴(𝑦) シーン表面の一部 𝑑𝜎(𝜔) 𝑑𝐴の法線 𝜃 𝑑𝜎(𝜔) = cos 𝜃 𝑟2 𝑑𝐴(𝑦) 𝜔 𝜔 = 𝑟 𝑦

Slide 35

Slide 35 text

確率密度の変換 • 微小立体角と微小面積の関係式より、立体角測度に関する確率密度と面積測度に関する確率 密度の関係式が得られる。以下の式により、確率密度の測度を変換することができる。 𝑃𝐴 (𝑦) = cos 𝜃 𝑟2 𝑃𝜎 (𝜔) 𝑃𝐴 𝑦 𝑑𝐴(𝑦) 𝑀 = 1, 𝑃𝜎 (𝜔)𝑑𝜎(𝜔) Ω = 1 𝑑𝜎(𝜔) = cos 𝜃 𝑟2 𝑑𝐴(𝑦) 𝑃𝜎 (𝜔)𝑑𝜎(𝜔) Ω = 𝑃𝜎 (𝜔) cos 𝜃 𝑟2 𝑑𝐴(𝑦) 𝑀 = 𝑃𝐴 𝑦 𝑑𝐴(𝑦) 𝑀 より （確率密度関数は全範囲で積分すると1） 導出

Slide 36

Slide 36 text

確率密度の変換（まとめ） 𝑥 𝜔 𝑦  位置𝑥から半球上に方向をサンプリングしてレイトレ―シング、その交点をサ ンプリングするということは、位置𝑥から見える全ての範囲（下図の赤い領 域）の中から一点をサンプリングするということと同じ。  前者は立体角測度に関する確率密度に基づいて𝑦をサンプリングしており、後 者は面積測度に関する確率密度に基づいて𝑦をサンプリングしている、という ことになる。  両者を相互に変換するのが先の式。

Slide 37

Slide 37 text

確率密度の変換（まとめ）  𝑷𝑨 (𝒚) = 𝒄𝒐𝒔 𝜽 𝒓𝟐 𝑷𝝈 (𝝎)  位置𝑥から半球上に方向をサンプリングするとき、 𝑥から近いシーン内の領域は 𝑟が小さく、その領域から𝑦がサンプリングされる確率（＝ 𝑥から発射したレイ がヒットする確率）も高くなる。逆に、 𝑥から遠いシーン内の領域は𝑟が大きく、 その領域から𝑦がサンプリングされる確率は低くなる。  また、対象の領域が傾いていればいるほど、レイがヒットする確率は下がる。 （＝cos 𝜃が小さくなる）  という関係を上の式は示している。  cos 𝜃 𝑟2 は、 位置𝑥からの𝑑𝐴(𝑦)の見かけの面積を計算している。 𝑥 𝜔 𝑦

Slide 38

Slide 38 text

具体例  𝑓 𝜔 𝑑𝜎(𝜔) Ω を解く際、変数変換して積分範囲や積分変数を変換すること を考える。  シーン内における微小立体角と微小面積の関係式𝑑𝜎 𝜔 = 𝑑𝜎(𝑥 → 𝑥′) = cos 𝜃 𝑟2 𝑑𝐴(𝑥′)より  𝑓 𝜔 𝑑𝜎(𝜔) Ω = 𝑓 𝑥 → 𝑥′ 𝑑𝜎(𝑥 → 𝑥′) Ω = 𝑓 𝑥 → 𝑥′ cos 𝜃 𝑟2 𝑑𝐴(𝑥′) M  ただしMは𝑥から一方向で到達可能な全てのシーンジオメトリ表面になる。  これをモンテカルロ積分するとなると、 M上の確率密度関数PA に基づいて サンプルを生成することになる  条件を満たしさえすれば任意の関数が与えられる – 面積測度に関する確率密度関数を使ったモンテカルロ積分に帰着することは多 いが、実際のレンダリングの文脈ではシーンジオメトリ上の確率密度関数を与 えて、それに基づいてサンプリングを行うことは少なく、結局立体角測に関す る確率密度関数によってサンプリングを行ってから確率密度関数を面積測度関 するものに変換することが多い。（例外：光源上のサンプリング、レンズ上の サンプリング）

Slide 39

Slide 39 text

1.4 ジオメトリファクタ

Slide 40

Slide 40 text

ジオメトリファクタ  シーン内の二頂点間の幾何的関係を示す項として、ジオメトリファクタ （以下、Ｇ項）と呼ばれるものがある。  今後、頻繁に出てくるので注意。 𝑥 𝑦 𝜃𝑥 𝜃𝑦 距離 𝑟 𝐺 𝑥 ↔ 𝑦 = cos 𝜃𝑥 cos 𝜃𝑦 𝑟2

Slide 41

Slide 41 text

1.5 マルチプルインポータンスサンプリング

Slide 42

Slide 42 text

マルチプルインポータンスサンプリング  ある積分をモンテカルロ積分で解くとする  普通のモンテカルロ積分は以下のようになる – このとき、サンプル𝑋𝑖 は確率密度関数𝑝(𝑥)に従ってサンプリングされている 𝐼 = 𝑓(𝑥)𝑑𝜇(𝑥) Ω 𝐼 = 1 𝑁 𝑓 𝑋𝑖 𝑝(𝑋𝑖 ) 𝑁 𝑖=1

Slide 43

Slide 43 text

マルチプルインポータンスサンプリング  モンテカルロ積分の効率（収束の速さ）はサンプリングのための確率密度 関数𝒑(𝒙)の形状によって決まる。 – 積分対象の関数𝑓(𝑥)に近ければ近いほど良い。 – 極論、 𝑓(𝑥)に比例していれば𝑓 𝑋𝑖 𝑝(𝑋𝑖) が定数になるため誤差は常に０になる。が、そもそも𝑓 𝑥 の形状が非常に複雑だからモンテカルロ積分を試みているため、これは不可能。 – それでも𝑓(𝑥)についての知識を利用してうまい確率密度関数を考えることは出来る。→ イ ンポータンスサンプリング（Importance Sampling） – 例えば、拡散面の上においてレンダリング方程式に出現するcos 𝜃に着目してcos 𝜃に比例し て次のレイの方向を決める、といったことが行われる。 𝐼 = 1 𝑁 𝑓 𝑋𝑖 𝑝(𝑋𝑖 ) 𝑁 𝑖=1

Slide 44

Slide 44 text

マルチプルインポータンスサンプリング  複数の確率密度関数を用意して、それぞれのモンテカルロ積分の結果を組 み合わせるという発想。  個々の確率密度関数は最適なものではないが、組み合わせれば良い結果が 得られえる → マルチプルインポータンスサンプリング（Multiple Importance Sampling） – 例えば、ある材質上で次のレイの方向をサンプリングするとき、(1)材質のBRDFに基づいて 確率密度関数を決めてサンプリング (2)光源の方向に基づいて確率密度関数を決めてサンプ リング (3)何も考えずランダムにサンプリング … といったように、複数のサンプリング戦略 が考えられる。このとき、唯一つの確率密度関数を上手に作成するのも良いが、そうではな く、別々に用意した確率密度関数を使って独立にサンプリングしてそれぞれのモンテカルロ 積分の結果を組み合わせる、といった戦略もあり得る。ということ。  以下のように定式化される。 – 𝑆は異なるサンプリング戦略（確率密度関数）の数。 – 𝑤𝑠 は各サンプリング戦略による各サンプルごとの重み。 𝐼 = 1 𝑁𝑠 𝑤𝑠 (𝑋𝑠,𝑖 ) 𝑓 𝑋𝑠,𝑖 𝑝𝑠 (𝑋𝑠,𝑖 ) 𝑁𝑠 𝑖=1 𝑆 𝑠=1

Slide 45

Slide 45 text

マルチプルインポータンスサンプリング  モンテカルロ積分がUnbiasedであるための重み𝒘𝒔 の条件 その１  モンテカルロ積分がUnbiasedであるための重み𝒘𝒔 の条件 その２ 𝑤𝑠 𝑥 ≠ 0 ⟹ 𝑝𝑠 𝑥 ≠ 0 あるいは 𝑝𝑠 𝑥 = 0 ⟹ 𝑤𝑠 𝑥 = 0 𝑓 𝑥 ≠ 0 ⟹ 𝑤𝑠 (𝑥) 𝑆 𝑠=1 = 1 以上は、任意のサンプル𝑥について、 𝑓(𝑥)が0でないならそのサンプルをサンプリ ングできる確率密度関数が少なくとも一つは存在しなければならず、さらにその ときの重みの和は1でなければならない、ということを意味する。

Slide 46

Slide 46 text

マルチプルインポータンスサンプリング  条件を満たす重みの例 – 全てのサンプリング戦略（確率密度関数）に対して等しい重みを与える。（た だし、任意のサンプルについて、全ての確率密度関数がサンプル出来なければ ならない） – あまり良くない。 𝑤𝑠 𝑥 = 1 𝑆

Slide 47

Slide 47 text

マルチプルインポータンスサンプリング  条件を満たす重みの例 – バランスヒューリスティック – 他の情報が無ければ最良の重み付け戦略となる。  これは、あるサンプル𝑥がサンプリング戦略𝑠によってサンプリングされたときの重 み。 – 𝑥が仮に他のサンプリング戦略によってサンプリングされたとして、その確率密 度の和に対する、 𝑠による確率密度の割合が上式。 𝑤𝑠 𝑥 = 𝑝𝑠 (𝑥) 𝑝𝑗 (𝑥) 𝑆 𝑗=1

Slide 48

Slide 48 text

具体例  二つのサンプリング戦略𝑝1 (𝑥)と𝑝2 (𝑥)があるとする。今、確率密度関数の 範囲をD = {𝑋, 𝑌}とすると、これらは以下の値をとるものとする。 – 𝑝1 𝑋 = 0.2, 𝑝1 𝑌 = 0.8 – 𝑝2 𝑋 = 0.6, 𝑝2 𝑌 = 0.4  今、 𝑝1 (𝑥)を使ってサンプリングした所、 𝑌が得られたとする。このときの 重み𝑤1 𝑌 は – ここで、 𝑌が仮に𝑝2 によってサンプリングされたとして、 𝑝2 (𝑌)を使っていることに注意。 𝑤1 𝑌 = 𝑝1 (𝑌) 𝑝1 𝑌 + 𝑝2 (𝑌) = 0.8 0.8 + 0.4 = 0.666 …

Slide 49

Slide 49 text

具体例  同様にして、全ての重みは以下のようになる。 – これは先ほどの条件を満たしている。  このことから、バランスヒューリスティックによる重み付けは、そのサンプルをサ ンプリングする確率密度が高い戦略に、より大きい重みが付けられるということが 分かる。  もし、一般のモンテカルロ積分のように𝑝1 のみ使ってサンプリングした場合、𝑋があ まりサンプリングされない。もし𝑓(𝑋)の値が大きい場合、これは大きな誤差・分散 を生む。また、𝑝2 のみの場合も逆のことが言える。そこで、二つの確率密度関数を 使い、それぞれの得意なサンプルについての重みを大きくすることで、ノイズの大 きくなるかもしれない結果の重みを小さくしつつ、ノイズが小さくなることが見込 まれる結果の重みを大きくしている。これがバランスヒューリスティックの考え方。 𝑤1 𝑋 = 0.25, 𝑤1 𝑌 = 0.666 … 𝑤2 𝑋 = 0.75, 𝑤2 𝑌 = 0.333 …

Slide 50

Slide 50 text

他の重みの例  条件を満たす重みの例 – パワーヒューリスティック – バランスヒューリスティックの一般化 – edubptではβ=2として使用。 𝑤𝑠 𝑥 = 𝑝𝑠 𝑥 𝛽 𝑝𝑗 𝑥 𝛽 𝑆 𝑗=1

Slide 51

Slide 51 text

2. レンダリングの定式化

Slide 52

Slide 52 text

レンダリングの定式化  レンダリングという問題を解くために、まずはどのような式でレン ダリングという問題を表現できるかをはっきりさせる。  そのために、カメラモデルと画素の受ける光について定式化する。

Slide 53

Slide 53 text

カメラモデル  イメージセンサに入射する光を物理的に計算し、それを各画素値に変換す ることで最終的な画像を得る。→ レンダリング  一般的なカメラの物理構造を模倣し、画像の計算モデルとして使う。 カメラモデル イメージセンサ レンズ レンダリング対象のシーン 光源

Slide 54

Slide 54 text

カメラモデル  イメージセンサ上の𝑗番目の画素の値について考える カメラモデル イメージセンサ レンズ 光源 画素𝑗 レンダリング対象のシーン

Slide 55

Slide 55 text

画素の受ける光  画素𝑗の最終的な値を𝐼𝑗 とする。  この値は、画素𝑗における放射束（単位時間あたりに画素を通過するフォト ン数・エネルギー）と画素のセンシティビティ（感応度）によって決まる とする。 – 今回は露光時間については考慮せず、単位時間あたりのエネルギーについて考える。 イメージセンサ 画素𝑗 𝐼𝑗 = 𝑊 𝑥𝐼 , 𝜔 𝐿 𝑥𝐼 , 𝜔 cos 𝜃 Ω 𝑑𝜎 𝜔 𝐼𝑗 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) センシティビティ 𝑥𝐼 における入射放射輝度×コサイン項 イラディアンスを面積分しているので得られるのは放射束に なる。（ただしセンシティビティが入る） 𝑥𝐼 𝑥𝐼 に𝜔方向から入射する放射輝度 𝜃 𝐿 𝑥𝐼 , 𝜔

Slide 56

Slide 56 text

光の波長  光は物理量としては波長ごとの放射輝度や放射照度として表現される。一 方、画像に実際に記録されるのはRGB値である。これは、RGBの三種類だけ で人間の認識できる色を表現できるからである。  人間の網膜には色を知覚する錐体細胞がある。この錐体細胞の種類は三つ あることが知られており、その三つそれぞれが波長ごとに異なった刺激を 受ける。つまり、人間の色の知覚は全てこの三種類の細胞の刺激の組み合 わせで表現できる。  そこで、三つの異なる波長の光をRGB原刺激として定め、これを様々な強 度で組み合わせることによって任意の色を表現することが行われている。 – たとえば、R刺激として波長700nm、G刺激として546.1nm、B刺激として435.8nmの光を考え、 この三種類の刺激を組み合わせることで任意の色を表現する（CIE 1931）。

Slide 57

Slide 57 text

光の波長  波長ごとの色の知覚と、対応するRGB原刺激値は以下のようになる。  以下のグラフは実験的に定めたものなので、実験した年や環境によっても 変化することに注意。 -1.00E+00 -5.00E-01 0.00E+00 5.00E-01 1.00E+00 1.50E+00 2.00E+00 2.50E+00 3.00E+00 3.50E+00 390 410 430 450 470 490 510 530 550 570 590 610 630 650 670 690 710 730 750 770 790 810 830 R G B 波長ごとのRGB刺激値（Stiles & Burch (1959) 10-deg, RGB CMF） 波長（nm） 刺激値

Slide 58

Slide 58 text

光の波長  人間が知覚できる色を表現するのはRGB値で十分であった。そこで、レン ダリングするときも全波長分の光の挙動を考慮するのではなく、RGB値に 対応する波長のみ考慮する。Color構造体が三要素なのも、そのため。  また、屈折率は波長ごとに変化するが、これも一種類とする。  この簡易化により、光の波長の連続的な変化に由来する現象は再現できな くなる。（光の分散等）

Slide 59

Slide 59 text

センシティビティ  RGBの三波長分のみ考慮してレンダリングするのだが、イメージセンサが 受ける光の量というのはあくまで物理量である。  人間の網膜上では、入射した光は最終的に「刺激」という心理的な量にな る。デジタルカメラのイメージセンサも、受け取る光は物理的実体だが、 最終的に画像に保存するときには何らかのスケーリングが行われRGBにつ いて相対的な量に変換される。この変換の係数がセンシティビティ。  今回のレンダラでは、センシティビティは画素の物理的大きさに反比例さ せている。このようにすることで、イメージセンサ全体の大きさを固定し て、解像度だけ変えても、各画素の値が変化しないようにしている。 – もしこうしないと解像度が変化すると各画素の物理的大きさが変化し、受け取る光の量も変 わり、画素値も変化するため不便。

Slide 60

Slide 60 text

画素の受ける光  画素𝑗内の点𝑥𝐼 において半球積分するわけだが、実際はレンズを通ってきた 光のみが画素𝑗に寄与する。  半球積分する代わりに、レンズ上の点𝑥0 から𝑥𝐼 への放射輝度𝐿 𝑥𝐼 ← 𝑥0 を 積分することで𝑥𝐼 におけるイラディアンスを計算、最終的な画素値を求め ることを考える。 – 𝑥0 をレンズ上を動く積分変数にする。  これは、積分変数を変数変換すればよい。 – 𝜔から𝑥0 に変換する。 イメージセンサ レンズ レンズモデルとして、厚さのない平面でレンズを近似する薄レンズモデルを使う。 𝑥𝐼 𝑥0 光軸 𝐿 𝑥𝐼 ← 𝑥0

Slide 61

Slide 61 text

画素の受ける光  微小立体角と微小面積の関係式𝑑𝜎(𝜔) = cos 𝜃 𝑟2 𝑑𝐴 𝑥0 を使うことで積分変数 と積分範囲を変換できる。 – 今回は、レンズをシーン中のほかのオブジェクトと同じ物理的実体とみなすため、レンズ上 の微小面積はシーン表面に対する測度と同じ面積測度を用いて得ることが出来ると考える。 イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐿 𝑥𝐼 ← 𝑥0 𝜃 𝐼𝑗 = 𝑊 𝑥𝐼 , 𝜔 𝐿 𝑥𝐼 , 𝜔 cos 𝜃 Ω 𝑑𝜎 𝜔 𝐼𝑗 𝑑𝐴𝐼𝑗 𝑥𝐼 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥𝐼 ← 𝑥0 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 𝜃 距離 𝑟

Slide 62

Slide 62 text

画素の受ける光  𝐿 𝑥𝐼 ← 𝑥0 を計算するために、シーン内の点𝑥1 から𝑥0 に到達し、レンズを通 過して𝑥𝐼 へと至るパスを計算する必要がある。 – このとき、 𝐿 𝑥𝐼 ← 𝑥0 = 𝐿 𝑥0 ← 𝑥1 となる。 – 𝒙𝑰 と𝒙𝟏 はそれぞれイメージセンサ上の点、シーン上の点で、別なので注意。  このような点𝑥1 は、𝑥𝐼 と𝑥0 が定まると自然に決まる。 – レンズの公式・結像公式 イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐿 𝑥𝐼 ← 𝑥0 𝜃 𝜃 シーン 𝑥1 求めたい点 𝐿 𝑥0 ← 𝑥1

Slide 63

Slide 63 text

レンズの公式  幾何光学における、レンズと焦点の関係式。 – イメージセンサとレンズの距離をB、レンズとオブジェクトプレーンの距離をA、レンズの 焦点距離をFとすると、以下の関係式が成り立つ。 – ここでオブジェクトプレーンとは、イメージセンサ上の一点から放射された光が再び一点で 収束する点が作る平面のことをいう。また、逆も成り立つ。 イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐵 𝐴 1 𝐴 + 1 𝐵 = 1 𝐹 オブジェクトプレーン 𝑥𝑉

Slide 64

Slide 64 text

許容錯乱円  イメージセンサとオブジェクトプレーンの位置関係によっては、オブジェ クトプレーン上で点であった物体がイメージプレーン上では円になる。 （錯乱円） – この円が、イメージセンサ上での画素（受光素子）よりも大きくなると、周囲の画素に影響 が波及し、いわゆるピンボケが起こる。 – この限界の円を許容錯乱円という。 イメージセンサ レンズ 光軸 錯乱円 オブジェクトプレーン

Slide 65

Slide 65 text

被写界深度  錯乱円の大きさが許容錯乱円以下になる（＝焦点が合う）、オブジェクト プレーンの範囲を被写界深度と呼ぶ。 イメージセンサ レンズ 光軸 被写界深度 オブジェクトプレーン

Slide 66

Slide 66 text

ボケ（Bokeh）  今回はレンズ形状を円形にしているため、錯乱円の形状も円形。  物理的なカメラは絞りと呼ばれる機構によって入射する光の量を制御して いる。  絞りの形状によって見かけのレンズ形状が変化するため、錯乱円の形状が 円形以外（六角形、五角形等）を取ることがある。 – 今回は行っていないが、レンズ上で𝑥0 をサンプリングする際、絞り形状の中からサンプリン グするようにすれば、任意の形状のボケを得ることが出来る。 今回はカメラモデルとして有限の大きさを持つレンズの存在 を考慮しているため、これらのボケ効果や被写界深度による ピンボケなどは全て再現される。

Slide 67

Slide 67 text

画素の受ける光  以上より、イメージセンサとレンズの距離、レンズの焦点距離が決まると レンズとオブジェクトプレーンの距離・位置も決まる。  よって、 𝑥𝐼 とレンズの中心点を結ぶ直線とオブジェクトプレーンの交点を 求めることが出来る。これを𝑥𝑉 とする。  𝑥0 から𝑥𝑉 に向かってレイトレーシングしてシーンとの交点を求めれば、そ れが𝑥1 となる。 イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐵 𝐴 オブジェクトプレーン 𝑥𝑉 シーン 𝑥1

Slide 68

Slide 68 text

画素の受ける光  以上より、画素𝑗の値は  𝐿 𝑥0 ← 𝑥1 はレンダリング方程式を使えば計算可能。 – パストレーシング イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐵 𝐴 オブジェクトプレーン 𝑥𝑉 シーン 𝑥1 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 𝐿 𝑥0 ← 𝑥1 𝐿 𝑥𝐼 ← 𝑥0

Slide 69

Slide 69 text

ビネッティング  以下は、白い光源だけをレンダリングした結果。画像周辺部分が暗くなっ ている。  これは、現実のカメラにおけるビネッティング（周辺光量低下）である。 特に、ナチュラルビネッティングと呼ぶ。  ひとつ前のスライドの積分中、 cos2 𝜃 𝑟2 について、𝑟 = 𝐵 cos 𝜃 で置き換えると 、 cos4 𝜃 𝐵2 となる。所謂コサイン四乗則である。  イメージセンサの端に行くほど、 𝜃の最大値も大きくなり、コサイン四乗 則による減衰も大きくなるため、画像周辺部が暗くなる。

Slide 70

Slide 70 text

3. パストレーシング

Slide 71

Slide 71 text

パストレーシング  前章で解くべき式が得られた。  双方向パストレについて考える前に、この式を普通のパストレーシングで 解いてみる。 – reference_pathtracing.hがこれに対応している。 – カメラモデル以外は、eduptのパストレーシングと大きく変わらない。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗

Slide 72

Slide 72 text

レンダリング方程式  イメージセンサの式とレンダリング方程式の二つがあれば問題を解くこと ができる。 – 特に、下のレンダリング方程式についてはeduptにおけるやり方と全く同様にし て解くことができる。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 𝐿 𝑥0 ← 𝑥1 = 𝐿𝑒 𝑥0 ← 𝑥1 + 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 cos 𝜃 𝑑𝜎(𝑥1 ← 𝑥2 ) Ω BRDF レンダリング方程式

Slide 73

Slide 73 text

対応するソースコード  camera.hのstruct Cameraに、カメラパラメータを設定。  render.hのrender_by_refernce_pathtracing()が画像の各画素値を計算する関数 で、𝐼𝑗 を解いている。 – 最初の点、 𝑥0 と𝑥𝐼 をCameraを使ってサンプリング。  render_by_refernce_pathtracing()の中で呼んでいるradiance_no_recursion がレ ンダリング方程式を解く関数で、パストレーシング部分。

Slide 74

Slide 74 text

render.h  render_by_reference_pathtracing() – この関数が画素値についてのモンテカルロ積分を行い、実際にレンダリングを 担当。

Slide 75

Slide 75 text

render.h  画素値についてのモンテカルロ積分 – レンズ上の点𝑥0 とイメージセンサ上の点𝑥𝐼 をサンプリング。 – 各点のサンプリング確率密度も得る。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 解いている式

Slide 76

Slide 76 text

render.h  画素値についてのモンテカルロ積分 – 続いて、𝐿 𝑥0 ← 𝑥1 を計算する。これはradiance_no_recursionが担当しており、 単純なパストレーシング。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 解いている式

Slide 77

Slide 77 text

render.h  画素値についてのモンテカルロ積分 – cos2 𝜃 𝑟2 を計算。 – 最後に𝑊 𝑥𝐼←𝑥0 𝐿 𝑥0←𝑥1 cos2 𝜃 𝑟2 𝑃𝐼𝑗 𝑥𝐼 𝑃𝐴(𝑥0) を計算して画素値に加算する（モンテカルロ積分）。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 解いている式

Slide 78

Slide 78 text

render.h  画素値についてのモンテカルロ積分 – 画像の出力時に、サンプリング数で割る。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 解いている式

Slide 79

Slide 79 text

camera.h • sample_points() • レンズ上の点𝒙𝟎 とイメージセンサ上の点𝒙𝑰 をサンプリング。 • まずイメージセンサ上で普通にサンプリング。 • その後、対応するオブジェクトプレーン上の位置を計算。 – これは前章で説明したとおり、レンズの公式を使って求まる。 • 最後に、レンズ上でサンプリング（円内にサンプリング点を生成）。

Slide 80

Slide 80 text

reference_pathtracing.h  radiance_no_recursion()によって𝐿 𝑥0 ← 𝑥1 を計算 – 再帰を使っていないだけでeduptのradiance関数とほぼ同じ。 – パストレーシングによってシーン内の光の経路を追跡する。

Slide 81

Slide 81 text

sampling.h  各種サンプリングアルゴリズムを提供 – インポータンスサンプリング用のcos 𝜃 𝜋 に従う半球サンプリング用関数等

Slide 82

Slide 82 text

specular.h  反射ベクトルの方向を得る関数や、屈折方向を得る関数を提供

Slide 83

Slide 83 text

4. 双方向パストレーシング・理論編

Slide 84

Slide 84 text

パストレーシング  カメラからレイを飛ばし、光源に当たるまで反射・追跡を繰り返す。これ によって、光源からカメラまでのパスを逆方向から得る。パスが得られれ ば、既に立式した積分に基づいて最終的な画像が計算できる。 𝑥0 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6

Slide 85

Slide 85 text

ライトトレーシング  光源からレイを飛ばし、カメラ（今回のカメラモデルではレンズ）に当た るまで反射・追跡を繰り返す。これによって、光源からカメラまでのパス を得る。パスが得られれば、既に立式した積分に基づいて最終的な画像が 計算できる。 𝑦6 𝑦5 𝑦4 𝑦3 𝑦2 𝑦1 𝑦0

Slide 86

Slide 86 text

双方向パストレーシング  カメラ側と光源側の両方からパスを追跡し、両者の結果を合わせることで 片側から追跡しただけではサンプリングしにくいパスを得る。 𝑥0 𝑥1 𝑥2 𝑥3 𝑦2 𝑦1 𝑦0

Slide 87

Slide 87 text

双方向パストレーシング  カメラ側と光源側の両方からパスを追跡し、両者の結果を合わせることで 片側から追跡しただけではサンプリングしにくいパスを得る。 𝑥0 𝑥1 𝑥2 𝑥3 𝑦2 𝑦1 𝑦0 新しいパス

Slide 88

Slide 88 text

4.1 双方向パストレが解く式

Slide 89

Slide 89 text

双方向パストレが解く式  画素𝑗の値について  であったが、双方向パストレーシングにおいてはこの式をさらに変形した ものを使う。（ライトトレーシングによって得られたパスを使う際、都合 が良くなるようにする）  具体的には、変数𝒙𝑰 の積分を変数𝒙𝟏 の積分に変数変換する。 – これによって、シーン上の点を変数とした積分になる。 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗

Slide 90

Slide 90 text

双方向パストレが解く式  イメージセンサ上の微小面積と、オブジェクトプレーン上の微小面積の関係式 – 相似関係が成り立つので、面積比は距離の二乗の比になる。 – よって、以下のようになる。 イメージセンサ レンズ 𝑥𝐼 光軸 𝐵 𝐴 オブジェクトプレーン 𝑥𝑉 𝑑𝐴𝐼𝑗 𝑥𝐼 𝑑𝑉 𝑥𝑉 𝑑𝐴𝐼𝑗 𝑥𝐼 : 𝑑𝑉 𝑥𝑉 = 𝐵2: 𝐴2 𝑑𝐴𝐼𝑗 𝑥𝐼 = 𝐵 𝐴 2 𝑑𝑉 𝑥𝑉

Slide 91

Slide 91 text

双方向パストレが解く式  まず、イメージセンサ上の微小面積と、オブジェクトプレーン上の微小面積の関 係式 𝑑𝐴𝐼𝑗 𝑥𝐼 = 𝐵 𝐴 2 𝑑𝑉 𝑥𝑉 を使い、変数𝑥𝐼 の積分を𝑥𝑉 の積分に変数変換する。 – このとき、 𝑊 𝑥𝐼 ← 𝑥0 = 𝑊 𝑥0 ← 𝑥𝑉 、𝐿 𝑥𝐼 ← 𝑥0 = 𝐿 𝑥0 ← 𝑥𝑉 という関係式を使う。  変換後の式 イメージセンサ レンズ 𝑥𝐼 𝑥0 光軸 𝐵 𝐴 オブジェクトプレーン 𝑥𝑉 𝑥1 𝐼𝑗 = 𝑊 𝑥0 ← 𝑥𝑉 𝐿 𝑥0 ← 𝑥𝑉 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝐵 𝐴 2 𝑑𝑉 (𝑥𝑉 ) 𝑉

Slide 92

Slide 92 text

双方向パストレが解く式  オブジェクトプレーン上の微小面積と、シーン上の微小面積の関係式 – 𝑥0 と𝑥𝑉 の距離を𝑟′′、 𝑥0 と𝑥1 の距離を𝑟′とする。 – 𝑑𝑉(𝑥𝑉 )をcos 𝜃′′で投影した面積と𝑑𝐴(𝑥1 )をcos 𝜃′で投影した面積の比が距離の二乗比になる。（相 似関係） レンズ 𝑥0 オブジェクトプレーン 𝑥𝑉 𝑥1 𝑑𝑉 𝑥𝑉 𝑑𝐴 𝑥1 𝜃′ 𝜃′′ 𝑟′′ 𝑟′ 𝑑𝑉 𝑥𝑉 cos 𝜃′′ : 𝑑𝐴 𝑥1 cos 𝜃′ = 𝑟′′2: 𝑟′2

Slide 93

Slide 93 text

双方向パストレが解く式  オブジェクトプレーン上の微小面積と、シーン上の微小面積の関係式 – 先の式より、  この関係式を使って変数𝑥𝑉 の積分を𝑥1 の積分に変数変換する。 – このとき、 𝑊 𝑥0 ← 𝑥𝑉 = 𝑊 𝑥0 ← 𝑥1 、𝐿 𝑥0 ← 𝑥𝑉 = 𝐿 𝑥0 ← 𝑥1 という関係式を使う。  この式をジオメトリファクタ（以下、Ｇ項）を使って書きなおす。 𝑑𝑉 𝑥𝑉 = 𝑟′′ 𝑟′ 2 cos 𝜃′ cos 𝜃′′ 𝑑𝐴(𝑥1 ) 𝐼𝑗 = 𝑊 𝑥0 ← 𝑥1 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝐵 𝐴 2 𝑟′′ 𝑟′ 2 cos 𝜃′ cos 𝜃′′ 𝑑𝐴(𝑥1 ) 𝑀 𝐺 𝑥0 ↔ 𝑥1 = cos 𝜃′ cos 𝜃′′ 𝑟′2

Slide 94

Slide 94 text

双方向パストレが解く式  すると以下のようになる。  𝑊′(𝑥0 ← 𝑥1 ) = 𝑊 𝑥0 ← 𝑥1 𝐵𝑟′′ cos 𝜃 𝐴𝑟 cos 𝜃′′ 2 とすると以下のようになる。 – シーンの面積に関する積分になった。 – この式を双方向パストレーシングで解く。 𝐼𝑗 = 𝑊 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐵𝑟′′ cos 𝜃 𝐴𝑟 cos 𝜃′′ 2 𝐿 𝑥0 ← 𝑥1 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴(𝑥1 ) 𝑀 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿 𝑥0 ← 𝑥1 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴(𝑥1 ) 𝑀

Slide 95

Slide 95 text

式変形のまとめ 𝐼𝑗 = 𝑊 𝑥𝐼 ← 𝑥0 𝐿 𝑥0 ← 𝑥1 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴𝐼𝑗 (𝑥𝐼 ) 𝐼𝑗 𝐼𝑗 = 𝑊 𝑥0 ← 𝑥𝑉 𝐿 𝑥0 ← 𝑥𝑉 cos2 𝜃 𝑟2 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝐵 𝐴 2 𝑑𝑉 (𝑥𝑉 ) 𝑉 𝐼𝑗 = 𝑊 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐵𝑟′′ cos 𝜃 𝐴𝑟 cos 𝜃′′ 2 𝐿 𝑥0 ← 𝑥1 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴(𝑥1 ) 𝑀 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿 𝑥0 ← 𝑥1 𝑑𝐴 𝑥0 𝐿𝑒𝑛𝑠 𝑑𝐴(𝑥1 ) 𝑀 𝑑𝐴𝐼𝑗 𝑥𝐼 = 𝐵 𝐴 2 𝑑𝑉 𝑥𝑉 𝑑𝑉 𝑥𝑉 = 𝑟′′ 𝑟′ 2 cos 𝜃′ cos 𝜃′′ 𝑑𝐴(𝑥1 ) 𝑊′(𝑥0 ← 𝑥1 ) = 𝑊 𝑥0 ← 𝑥1 𝐵𝑟′′ cos 𝜃 𝐴𝑟 cos 𝜃′′ 2 元の式 最終的な式

Slide 96

Slide 96 text

4.2 レンダリング方程式

Slide 97

Slide 97 text

レンダリング方程式  𝐿 𝑥0 ← 𝑥1 はレンダリング方程式を使って再帰的に定義される。 𝐿 𝑥0 ← 𝑥1 = 𝐿𝑒 𝑥0 ← 𝑥1 + 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 cos 𝜃 𝑑𝜎(𝑥1 ← 𝑥2 ) Ω BRDF 𝑥1 𝑥0 𝑥2 𝜃 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝜃′

Slide 98

Slide 98 text

レンダリング方程式  これは立体角に関する半球積分だが、シーン上の点を変数とする、面積に 関する積分に変換する。（積分範囲はシーン表面𝑀） – 𝐼𝑗 の積分を変数変換してシーン上の点を変数とする面積分にしたのと同じ理由。 – 例によって微小立体角と微小面積の関係式𝑑𝜎(𝑥1 ← 𝑥2) = cos 𝜃′ 𝑟2 𝑑𝐴 𝑥2 を使う。 – さらに𝐺 𝑥1 ↔ 𝑥2 = cos 𝜃 cos 𝜃′ 𝑟2 も使う。 𝐿 𝑥0 ← 𝑥1 = 𝐿𝑒 𝑥0 ← 𝑥1 + 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 cos 𝜃 𝑑𝜎(𝑥1 ← 𝑥2 ) Ω = 𝐿𝑒 𝑥0 ← 𝑥1 + 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝐺 𝑥1 ↔ 𝑥2 𝑑𝐴(𝑥2 ) 𝑀 𝑥1 𝑥0 𝑥2 𝜃 𝐿 𝑥1 ← 𝑥2 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝜃′

Slide 99

Slide 99 text

レンダリング方程式  すると、解く対象の積分は以下のように展開される。 𝑥1 𝑥0 𝐿𝑒 𝑥2 ← 𝑥3 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝑓𝑟 𝑥1 ← 𝑥2 ← 𝑥3 𝑥2 𝑥3 𝐺 𝑥1 ↔ 𝑥2 𝐺 𝑥2 ↔ 𝑥3 𝑊′ 𝑥0 ↔ 𝑥1 レンズ 𝐺 𝑥0 ↔ 𝑥1 光源

Slide 100

Slide 100 text

レンダリング方程式  つまり、光源からレンズまでのパスの頂点の数の分だけ、独立した重積分 が得られ、その和が最終的な画素値𝐼𝑗 になる、ということ。 𝑥1 𝑥0 𝐿𝑒 𝑥2 ← 𝑥3 𝑓𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝑓𝑟 𝑥1 ← 𝑥2 ← 𝑥3 𝑥2 𝑥3 𝐺 𝑥1 ↔ 𝑥2 𝐺 𝑥2 ↔ 𝑥3 𝑊′ 𝑥0 ↔ 𝑥1 レンズ 𝐺 𝑥0 ↔ 𝑥1 光源

Slide 101

Slide 101 text

モンテカルロ積分  そのような重積分一つ一つを、モンテカルロ積分で解くことが出来る。 – シーン上の頂点によって構成されるパスが、モンテカルロ積分に対するサンプルとして必要 になる。 – そのようなパスは、例えばパストレーシングなどで得ることが出来る。 – 双方向パストレーシングにおいては、カメラ側からのトレーシングと光源側からのトレーシ ングによって得られた頂点集合を組み合わせることで様々な頂点数の複数のパスを生成し、 それらのパスを各モンテカルロ積分の計算に使用する。 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿𝑒 (𝑥0 ← 𝑥1 ) 𝑃𝐴 𝑥0 𝑃𝐴 (𝑥1 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝐿𝑒 (𝑥1 ← 𝑥2 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 ) + ⋯ 𝐼𝑗 = 頂点数２の積分＋頂点数３の積分＋頂点数４の積分＋ …

Slide 102

Slide 102 text

5. 双方向パストレーシング・実装編その１

Slide 103

Slide 103 text

双方向パストレーシング・実装編  カメラ側からパストレーシングしてサンプリングされた頂点集合と、光源 側からライトトレーシングしてサンプリングされた頂点集合をマージして 一つのパスにする。 – そのやり方について、実装を並べつつ解説。 𝑥0 𝑥1 𝑥2 𝑥3 𝑦2 𝑦1 𝑦0 新しいパス

Slide 104

Slide 104 text

最終的に必要なもの  最終的に以下のようなモンテカルロ積分の和を計算しなければならない  そこで、各頂点ごとの面積測度に関する確率密度（分母）とモンテカルロ スループット（分子）を求める必要がある。 – カメラ側からのパス、光源側からのパス、両方についてスループットと確率密 度を求める。 – カメラ側のパスと光源側のパスをつなげるときは、それぞれの確率密度とス ループットを乗算していくことになる。  この結果が、上の式のモンテカルロ積分の和を計算するための一つ一つの 項の推定値になる。（上の式は頂点数ごとに別々のモンテカルロ積分の和 になっており、その各項） 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿𝑒 (𝑥0 ← 𝑥1 ) 𝑃𝐴 𝑥0 𝑃𝐴 (𝑥1 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝐿𝑒 (𝑥1 ← 𝑥2 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 ) + ⋯

Slide 105

Slide 105 text

最終的に必要なもの 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿𝑒 (𝑥0 ← 𝑥1 ) 𝑃𝐴 𝑥0 𝑃𝐴 (𝑥1 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝐿𝑒 (𝑥1 ← 𝑥2 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝑓𝑟 (𝑥1 ← 𝑥2 ← 𝑥3 )𝐺 𝑥2 ↔ 𝑥3 𝐿𝑒 (𝑥2 ← 𝑥3 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 𝑥2 𝑃𝐴 (𝑥3 ) + ⋯ パストレーシングによってサンプリングされた パス𝑥0 𝑥1 のスループットと確率密度 ライトトレーシングによってサンプリングされた パス𝑥2 𝑥3 のスループットと確率密度 頂点𝑥1 と𝑥2 を新しくつないだので、改めてG項と BRDFを計算  以上のような計算を行うために、前半部分をパストレーシング、後 半部分をライトトレーシングによって得るわけである。

Slide 106

Slide 106 text

最終的に必要なもの 𝑥0 𝑥1 𝑥2 𝑥3 𝑦2 𝑦1 𝑦0 上の色を付けた項と対応するパス 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿𝑒 (𝑥0 ← 𝑥1 ) 𝑃𝐴 𝑥0 𝑃𝐴 (𝑥1 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝐿𝑒 (𝑥1 ← 𝑥2 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝑓𝑟 (𝑥1 ← 𝑥2 ← 𝑥3 )𝐺 𝑥2 ↔ 𝑥3 𝐿𝑒 (𝑥2 ← 𝑥3 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 𝑥2 𝑃𝐴 (𝑥3 ) + ⋯

Slide 107

Slide 107 text

5.1 パストレーシングによるパスのサンプリング

Slide 108

Slide 108 text

pathtracing.h  パストレーシングによる頂点のサンプリングを行う – generate_vertices_by_pathtracing ()

Slide 109

Slide 109 text

パストレーシングによるサンプリングの例 𝑥0  まず、レンズ上𝑥0 とイメージセンサ上𝑥𝐼 にサンプルを生成する。この時、 各頂点のサンプリング確率密度も一緒に求める。 – これが𝑃𝐴 𝑥0 になる。ソースコード中ではP_lensである。  シーンをトレースし、次々に頂点をサンプリングしていき、各頂点をリス トに追加していく。 – 頂点𝑥𝑁 を追加するとき、 𝑁までの確率密度の総計（ 𝑃𝐴 𝑥𝑖 𝑁 𝑖=0 ）とモンテカルロスループッ トの総計（G項とBRDFの積）が必要になるので、順次計算しておく。 – それぞれソースコード中ではtotal_pdf_AとMC_throughputという変数に格納しておく。 𝑥𝐼 𝑃𝐴 (𝑥0 ) レンズ 𝑃𝐼𝑗 (𝑥𝐼 ) イメージセンサ

Slide 110

Slide 110 text

パストレーシングによるサンプリングの例 𝑥1  既に見たように、 𝑥0 と𝑥𝐼 が決まると、シーンをレイトレーシングすること で最初のシーン上の頂点𝑥1 をサンプリングできる。 – レンズ上の点から、対応するオブジェクトプレーン上の点の方向にレイを飛ばす。  𝑃𝐴 𝑥1 を計算する必要があるが、これには以下の関係式を使う。 𝑑𝐴𝐼𝑗 𝑥𝐼 = 𝐵 𝐴 2 𝑑𝑉 𝑥𝑉 𝑑𝑉 𝑥𝑉 = 𝑟′′ 𝑟′ 2 cos 𝜃′ cos 𝜃′′ 𝑑𝐴(𝑥1 ) 𝑥0 𝑥𝐼 𝑃𝐴 (𝑥0 ) レンズ 𝑃𝐼𝑗 (𝑥𝐼 ) イメージセンサ

Slide 111

Slide 111 text

パストレーシングによるサンプリングの例  𝑑𝐴 𝑎 = 𝐶𝑑𝐵 𝑏 ⟺ 𝐶𝑃𝐴 𝑎 = 𝑃𝐵 𝑏 という関係より、以下の式が成り立つ。 この式を使ってイメージセンサ上の確率密度𝑃𝐼𝑗 (𝑥𝐼 )をシーンの面積測度に 関する確率密度𝑃𝐴 𝑥1 に変換する。 レンズ 𝑥0 オブジェクトプレーン 𝑥𝑉 𝑥1 𝑑𝑉 𝑥𝑉 𝑑𝐴 𝑥1 𝜃′ 𝜃′′ 𝑟′′ 𝑟′ 𝑃𝐴 𝑥1 = 𝐵 𝐴 2 𝑟′′ 𝑟′ 2 cos 𝜃′ cos 𝜃′′ 𝑃𝐼𝑗 (𝑥𝐼 )

Slide 112

Slide 112 text

パストレーシングによるサンプリングの例 𝑥1  camera.P_image_to_PA_x1()が実際の変換式を実装したものになっている。 𝑃𝐴 (𝑥1 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 113

Slide 113 text

パストレーシングによるサンプリングの例 𝑥1  なお、各頂点ごとに生成するかしないかをロシアンルーレットで判定する ことで無限に頂点が生成されるのを防いでいる。  最終的な頂点のサンプリング確率密度𝑃𝐴 (𝑥𝑖 )はロシアンルーレットの確率を 含んだものになる。 – ロシアンルーレットの確率は頂点ごとに独立に計算する。 𝑃𝐴 (𝑥1 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 114

Slide 114 text

パストレーシングによるサンプリングの例 𝑥1  モンテカルロスループットは通常はG項とBRDFを乗算するだけだが、𝑥1 に 関しては𝑊′ 𝑥0 ← 𝑥1 の分を余分に考慮する必要がある。 – センサのセンシティビティや、レンズに光が入射することによる幾何的な係数 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 115

Slide 115 text

パストレーシングによるサンプリングの例 𝑥1  その後、G項を計算し、新しい頂点を頂点リストに追加する。 – その頂点までの確率密度の総計と、モンテカルロスループットの総計も追加する。  ただし、新しい頂点が光源上にあった場合、そこで追跡を打ち切る。 – 今回光源は全て反射率0と仮定しており、これ以上の追跡は無駄であるため。もし反射率0と いう仮定を入れないなら更なる追跡が必要になる。 – この段階のトータルのコントリビューション（モンテカルロスループット/パスの確率密 度）も計算しておく。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 116

Slide 116 text

パストレーシングによるサンプリングの例 𝑥1  次の頂点をサンプリングするために、今の頂点を始点としてレイトレーシン グを行う。これは立体角測度に関する確率密度によるサンプリングとなる。  使用する確率密度は現在の頂点の材質によって決める。 – edupt同様、今回は完全拡散面、完全鏡面、完全ガラス面の三種類を考える。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 117

Slide 117 text

完全拡散面の場合 𝑥1  普通に半球上から次の方向をサンプリングする。 cos 𝜃 を考慮しているため、 確率密度はcos 𝜃 𝜋 となる。この確率密度は次の頂点𝑥2 の分になる。  さらに、スループットにBRDFの分を乗算する。完全拡散面なので、 𝜌 𝜋 とな る。（ 𝜌は反射率）このスループットは次の頂点𝑥2 の分になる。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑃𝜎 (𝑥1 → 𝑥2 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 118

Slide 118 text

完全鏡面の場合 𝑥1  完全鏡面の場合、次のレイの方向は一方向に決まる。  ただし、その場合でも立体角測度に関する確率密度関数とBRDFを求める必要 がある。 – モンテカルロ積分に投入する必要があるため。  そこで、Diracのδ関数を導入する。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑠 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 119

Slide 119 text

Diracのδ関数  完全鏡面のBRDFとして望ましい性質。 – 光が完全にある一方向のみに反射されること。つまり以下のような性質。  そのようなBRDF、𝑓𝑠 としてDiracの𝛿関数を採用する。 – この関数を使うことで、完全鏡面やガラス面も普通の完全拡散面に近い取扱い を行える。  レンダリングの文脈では以下のような関数𝛿を考える。 – 𝑥 ≠ 0 ⟹ 𝛿 𝑥 = 0 – 𝛿 𝑥 𝑑𝑥 ℝ = 1 𝐿 𝑥, 𝜔𝑜 = 𝐿 𝑥, 𝜔𝑖 𝑓𝑠 𝜔𝑖 → 𝜔𝑜 cos 𝜃 𝑑𝜎 𝜔𝑖 𝑆2 = 𝐿 𝑥, 𝜔𝑟 𝜔𝑜 𝜔𝑟 完全鏡面 𝜃 𝜃 𝐿 𝑥, 𝜔𝑜 𝐿 𝑥, 𝜔𝑟

Slide 120

Slide 120 text

Diracのδ関数  𝛿関数を使って𝑓𝑠 を以下のように書くことが出来る。  これは確かに、以下の式を満たす。 𝑓𝑠 𝜔𝑖 → 𝜔𝑜 = 𝛿𝜎 𝜔𝑖 − 𝜔𝑟 cos 𝜃 𝐿 𝑥, 𝜔𝑜 = 𝐿 𝑥, 𝜔𝑖 𝑓𝑠 𝜔𝑖 → 𝜔𝑜 cos 𝜃 𝑑𝜎 𝜔𝑖 𝑆2 = 𝐿 𝑥, 𝜔𝑟 𝜔𝑜 𝜔𝑟 完全鏡面 𝜃 𝜃 𝐿 𝑥, 𝜔𝑜 𝐿 𝑥, 𝜔𝑟

Slide 121

Slide 121 text

Diracのδ関数  この時、立体角測度に関する確率密度関数𝑃𝜎 (𝜔𝑖 )は以下のようになる。 𝑃𝜎 (𝜔𝑖 ) = 𝛿𝜎 𝜔𝑖 − 𝜔𝑟

Slide 122

Slide 122 text

完全鏡面の場合 𝑥1  以上より、now_sampled_pdf_omegaは𝛿𝜎 𝜔𝑖 − 𝜔𝑟 となり、スループットにはBRDFである 𝛿𝜎 𝜔𝑖−𝜔𝑟 cos 𝜃 を乗算すればよいということになる。しかし、𝛿𝜎 を数値的に表現することはできな い。  記号的に𝛿𝜎 の存在を記録してもよいが、最終的なモンテカルロ積分においてBRDFが分子に、 確率密度に分母が来て、一対一対応がとれることから、𝛿𝜎 は最終的にお互い打ち消しあうと いう事実を利用する。  つまり、𝛿𝜎 を抜いた分をnow_sampled_pdf_omegaとすればよいし、モンテカルロスループッ トも抜いた分で計算する。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑠 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 123

Slide 123 text

完全ガラス面の場合 𝑥1  基本的に完全鏡面の場合と同じになる。全反射の場合は全く同じだし、反射＋ 屈折が起こる際も、ロシアンルーレットで反射が選ばれればやはりまったく同 じになる。 – もちろん、ロシアンルーレットの確率をtotal_pdf_Aに乗算する必要がある。  屈折が起こる場合も、大きくは変わらないが、屈折前と屈折後で立体角が変化 することによる放射輝度値の変化（nnt2）をスループットの計算に入れる。 𝑊′ 𝑥0 ← 𝑥1 𝑃𝐴 (𝑥1 ) 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑠 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 124

Slide 124 text

パストレーシングによるサンプリングの例 𝑥1  次のループの最初で、新しい頂点がシーン上からサンプリングされる。 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥2 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝐺 𝑥0 ↔ 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 125

Slide 125 text

パストレーシングによるサンプリングの例 𝑥1  最初の頂点𝑥1 以外は、面積測度に関する確率密度𝑃𝐴 (𝑥𝑖 )を、その頂点をサ ンプリングするための立体角測度に関する確率密度𝑃𝜎 (𝑥𝑖−1 → 𝑥𝑖 )から変換 して求める。 – 変換には、関係式𝑃𝐴 (𝑥𝑖 ) = cos 𝜃 𝑟2 𝑃𝜎 (𝑥𝑖−1 → 𝑥𝑖)を使う。 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥2 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑃𝐴 (𝑥2 ) 𝐺 𝑥0 ↔ 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 126

Slide 126 text

パストレーシングによるサンプリングの例 𝑥1  以上の処理を（ロシアンルーレットで打ち切られるか光源にヒットするま で）繰り返すことで頂点を生成していく。  同時に、モンテカルロ積分に使うためのスループットと確率密度も計算し ていく。 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥2 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑃𝐴 (𝑥2 ) 𝑥3 𝑥4 𝑓𝑟 (𝑥1 ← 𝑥2 ← 𝑥3 ) 𝑓𝑟 (𝑥2 ← 𝑥3 ← 𝑥4 ) 𝑃𝐴 (𝑥3 ) 𝑃𝐴 (𝑥4 ) 𝐺 𝑥1 ↔ 𝑥2 𝐺 𝑥0 ↔ 𝑥1 𝐺 𝑥2 ↔ 𝑥3 𝐺 𝑥3 ↔ 𝑥4 𝑥0 𝑃𝐴 (𝑥0 ) レンズ

Slide 127

Slide 127 text

パストレーシングによるサンプリングの例 𝑥1  このようにして生成されたパスはモンテカルロ積分に使うことが出来る。 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥2 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑃𝐴 (𝑥2 ) 𝑥3 𝑥4 𝑓𝑟 (𝑥1 ← 𝑥2 ← 𝑥3 ) 𝑓𝑟 (𝑥2 ← 𝑥3 ← 𝑥4 ) 𝑃𝐴 (𝑥3 ) 𝑃𝐴 (𝑥4 ) 𝐺 𝑥1 ↔ 𝑥2 𝐺 𝑥0 ↔ 𝑥1 𝐺 𝑥2 ↔ 𝑥3 𝐺 𝑥3 ↔ 𝑥4 𝑥0 𝑃𝐴 (𝑥0 ) レンズ 𝐼𝑗 = 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝐿𝑒 (𝑥0 ← 𝑥1 ) 𝑃𝐴 𝑥0 𝑃𝐴 (𝑥1 ) + 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 )𝐺 𝑥1 ↔ 𝑥2 𝐿𝑒 (𝑥1 ← 𝑥2 ) 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 ) + ⋯

Slide 128

Slide 128 text

5.2 ライトトレーシングによるパスのサンプリング

Slide 129

Slide 129 text

lighttracing.h  ライトトレーシングによる頂点のサンプリングを行う – generate_vertices_by_lighttracing()

Slide 130

Slide 130 text

ライトトレーシングによるサンプリングの例  まず、光源上𝑦0 にサンプルを生成する。頂点のサンプリング確率密度も求める。 – これが𝑃𝐴 𝑦0 になる。ソースコード中ではpdf_A_on_lightである。 – 半径Ｒの球上から一様にサンプリングしているため、 𝑃𝐴 𝑦0 = 1 4𝜋𝑅2 になる。  パストレーシングの場合と同様に、シーンをトレースし、次々に頂点をサンプリン グしていく。 – 頂点𝑦𝑁 を追加するとき、 𝑁までの確率密度の総計（ 𝑃𝐴 𝑦𝑖 𝑁 𝑖=0 ）とモンテカルロスループットの総計（G項 とBRDFの積）が必要になるので、順次計算しておく。 – それぞれソースコード中ではtotal_pdf_AとMC_throughputという変数に格納しておく。 – この辺もパストレーシングの場合と同様。  なお、完全拡散光源なので𝐿𝑒 𝑦1 ← 𝑦0 は次の頂点𝑦1 をサンプリングしなくても決ま る。 光 源 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0

Slide 131

Slide 131 text

ライトトレーシングによるサンプリングの例  光源の次の頂点を生成するために、レイの方向をサンプリングする。これ は、完全拡散面上における方向のサンプリングと同じようにして行う。 光 源 𝑃𝜎 (𝑦1 ← 𝑦0 ) 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0

Slide 132

Slide 132 text

ライトトレーシングによるサンプリングの例  レイの方向に向かってレイトレーシングして次の頂点を得る。  この時、レンズとの交差判定も行う。 – レンズにヒットした場合、そこで追跡を打ち切る。 – 幾何的な計算により、レンズ上の点と、対応するイメージセンサ上の 点が求まる。 光 源 𝑥0 𝑥𝐼 レンズ イメージセンサ 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0

Slide 133

Slide 133 text

ライトトレーシングによるサンプリングの例  レンズ上の点𝑥0 の面積測度に関する確率密度を計算する。これは、𝑃𝜎 (𝑦1 ← 𝑦0 ) の測度を面積測度に変換することで行う。 – パストレーシングの時と同様にして測度を変換する。  また、G項も計算する。  さらに、イメージセンサのセンシティビティ諸々を含めた項、W’も計算する。 – この段階のトータルのコントリビューション（モンテカルロスループット/パスの確率密度）も 計算しておく。 光 源 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 レンズ イメージセンサ 𝑥0 𝑥𝐼 𝑃𝐴 𝑥0 𝐺 𝑥0 ↔ 𝑦0 𝑊′ 𝑥0 ← 𝑦0

Slide 134

Slide 134 text

ライトトレーシングによるサンプリングの例  飛ばしたレイがレンズとヒットしなかった場合、シーンとの交点を新しい 頂点とする。  ロシアンルーレットも行う。  さらに、新しい頂点の面積測度に関する確率密度𝑃𝐴 (𝑦1 )を、立体角測度に 関する確率密度𝑃𝜎 (𝑦1 ← 𝑦0 )から求める。  G項も求め、スループットに乗算する。 光 源 𝑦1 𝑃𝐴 𝑦1 𝐺 𝑦1 ↔ 𝑦0 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0

Slide 135

Slide 135 text

ライトトレーシングによるサンプリングの例  現在の頂点の材質ごとに、次の頂点をサンプリングするための方向のサン プリング方法も変わる。  これは、パストレーシングの時と基本的に同じ処理になる。 – パストレーシングのときと同様、スループットにBRDFを乗算していくことになる。 – ただし、nnt2については光の進行方向の違いにより、若干処理が変わる。 光 源 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑦1 𝑃𝐴 𝑦1 𝐺 𝑦1 ↔ 𝑦0 𝑓𝑟 (𝑦2 ← 𝑦1 ← 𝑦0 )

Slide 136

Slide 136 text

ライトトレーシングによるサンプリングの例  以上の処理を（ロシアンルーレットで打ち切られるかレンズにヒットする まで）繰り返すことで頂点を生成していく。  同時に、モンテカルロ積分に使うためのスループットと確率密度も計算し ていく。最初の部分以外とレンズ判定以外はパストレーシングとほぼ同じ になる。 光 源 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑦1 𝑃𝐴 𝑦1 𝐺 𝑦1 ↔ 𝑦0 𝑦2 𝑓𝑟 (𝑦2 ← 𝑦1 ← 𝑦0 ) 𝑓𝑟 (𝑦3 ← 𝑦2 ← 𝑦1 ) 𝑦3 𝑃𝐴 𝑦2 𝑃𝐴 𝑦3 𝐺 𝑦2 ↔ 𝑦1 𝐺 𝑦3 ↔ 𝑦2

Slide 137

Slide 137 text

まとめ  解きたいモンテカルロ積分について、 – 分子（モンテカルロスループット） • 光源における放射輝度×それ以降の頂点間のBRDFの積×それ以降の頂点間のG項の積× センサのセンシティビティ等の項（W’） – 分母（確率密度） • 各頂点をサンプリングするための面積測度に関する確率密度の積  が必要なので、パストレーシング or ライトトレーシングによってこれを頂 点ごとに逐次求めていく。 – この結果を組み合わせて双方向パストレーシングを実行する。

Slide 138

Slide 138 text

6. 双方向パストレーシング・実装編その２

Slide 139

Slide 139 text

パスの接続の例  パストレーシングによってサンプリングされたパス𝑥0 𝑥1 𝑥2 とライトトレー シングによってサンプリングされたパス𝑦1 𝑦0 を接続すると、以下のような 五頂点分のモンテカルロ積分を解くことに相当する。 – total_pdf_Aや、MC_throughputは、以下を計算するためのものだった。 – 下の式の緑の部分の分母と分子が、generate_vertices_by_pathtracing ()によって得られる頂点 リストの要素であるVertex構造体のメンバ、total_pdf_A とthroughputに対応する。 • この場合は、𝑥2 の分までなので、3番目の要素のメンバになる。 – 同様にオレンジの部分はgenerate_vertices_by_lighttracing ()によって得られる。 – 黒い部分は頂点の接続によって新たに計算される。 𝑥1 𝑃𝐴 (𝑥1 ) 𝑊′ 𝑥0 ← 𝑥1 𝑥2 𝑓 𝑟 (𝑥0 ← 𝑥1 ← 𝑥2 ) 𝑃𝐴 (𝑥2 ) 𝐺 𝑥1 ↔ 𝑥2 𝐺 𝑥0 ↔ 𝑥1 𝑥0 𝑃𝐴 (𝑥0 ) レンズ 光 源 𝑦0 𝑃𝐴 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑦1 𝑃𝐴 𝑦1 𝐺 𝑦1 ↔ 𝑦0 𝑓𝑟 (𝑥1 ← 𝑥2 ← 𝑦1 ) 𝑓𝑟 (𝑥2 ← 𝑦1 ← 𝑦0 ) 𝐺 𝑥2 ↔ 𝑦1 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓 𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝐺 𝑥1 ↔ 𝑥2 𝑓 𝑟 𝑥1 ← 𝑥2 ← 𝑦1 𝐺 𝑥2 ↔ 𝑦1 𝑓 𝑟 𝑥2 ← 𝑦1 ← 𝑦0 𝐺 𝑦1 ↔ 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 )𝑃𝐴 (𝑦0 )𝑃𝐴 (𝑦1 )

Slide 140

Slide 140 text

パスの接続とモンテカルロ積分  以下、具体的な例。  頂点数が3のカメラ側のパスと頂点数が2のライト側のパスを接続した場合、 新たに頂点数5のパスが得られる。  頂点数が2のカメラ側のパスと頂点数が3のライト側のパスを接続した場合 も、新たに頂点数5のパスが得られる。  この二つのパスは、頂点数5の積分に対するモンテカルロ積分を異なるサン プリング戦略でサンプリングして得られたサンプル（パス）である、と考 える。  よって、それぞれのパスによるモンテカルロ積分を独立に行い、その結果 をマルチプルインポータンスサンプリングによって重み付けして統合する 必要がある。 – 同じ積分に対するモンテカルロ積分を、異なるサンプリング戦略で独立に解いた場合、各々 の結果をマルチプルインポータンスサンプリングによって統合するのだった。 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥𝑖 ↔ 𝑥𝑖+1 3 𝑖=0 𝑓𝑟 𝑥𝑖 ← 𝑥𝑖+1 ← 𝑥𝑖+2 2 𝑖=0 𝐿𝑒 (𝑥3 ← 𝑥4 )𝑑𝐴(𝑥0 ) 𝐿𝑒𝑛𝑠 𝑑𝐴(𝑥1 ) 𝑀 𝑑𝐴(𝑥2 ) 𝑀 𝑑𝐴(𝑥3 ) 𝑀 𝑑𝐴(𝑥4 ) 𝑀 頂点数5の積分

Slide 141

Slide 141 text

光 源 光 源 光 源 光 源 光 源 光 源 同じ頂点数の同じ積分を解いているが、頂点の組み合わせの数だけ別のモン テカルロ積分が存在し、それらの結果をMISによって統合する。 カメラ側のパス：頂点数0 光源側のパス：頂点数5 カメラ側のパス：頂点数1 光源側のパス：頂点数4 カメラ側のパス：頂点数2 光源側のパス：頂点数3 カメラ側のパス：頂点数3 光源側のパス：頂点数2 カメラ側のパス：頂点数4 光源側のパス：頂点数1 カメラ側のパス：頂点数5 光源側のパス：頂点数0

Slide 142

Slide 142 text

bidirectional_pathtracing.h  頂点同士を接続し、最終的な画像を得るための実装を説明していく。

Slide 143

Slide 143 text

双方向パストレーシング  bidirectional_pathtracing()によって、カメラパラメータ、現在のピクセル座 標（𝐼𝑗 の𝑗に対応する 𝑥, 𝑦 ）を受け取り、双方向パストレーシングによって 画素値を計算する。  まず、パストレーシングとライトトレーシングによって、カメラ側、光源 側より頂点をサンプリングしてパスを生成する。

Slide 144

Slide 144 text

双方向パストレーシング  パストレーシングによるパスが直接光源にヒットしていた場合、その結果 をそのまま画素値の計算に使う。（この時、ライトトレーシングによる頂 点数は0）  また、ライトトレーシングによるパスが直接レンズにヒットしていた場合、 その結果をそのまま対応する画素値の計算に使う。（この時、パストレー シングによる頂点数は0）  既に説明したように、それぞれのパスによるモンテカルロ積分に対するMIS の重みを計算する。後で詳しく解説。

Slide 145

Slide 145 text

頂点間の接続  以後、カメラ側のパスの一部を使って作られるパスをeyeサブパス、光源側 のパスの一部を使って作られるパスをlightサブパスと呼ぶことにする。  eyeサブパスとlightサブパスを接続するのが今後の処理になる。全ての接続の 組み合わせを試す。 – これによってさまざまな頂点数のパスが得られ、さまざまな頂点数の積分に対するモンテカ ルロ積分の結果を得られる。  例えば、カメラ側のパスの頂点数が7、光源側のパスの頂点数が5だとすると、 eyeサブパスは頂点数1から7の7種類、lightサブパスは頂点数1から5の5種類が 考えられる。この時、接続の組み合わせは7×5=35種類となる。 – ひとつ前のスライドの処理は、eyeサブパスやlightサブパスの頂点数が0の場合の処理であった。

Slide 146

Slide 146 text

頂点間の接続  eyeサブパスの頂点数とlightサブパスの頂点数が決まったら、いよいよサブ パスを接続する。  まず、接続したパスのトータルのサンプリング確率密度を計算する。これ は、各サブパスのサンプリング確率密度を乗算すればよい。 – 接続したパスのサンプリング確率密度とは、パスを構成する頂点の面積測度に関する確率密 度の積であった。そして、各サブパスのサンプリング確率密度は、サブパスを構成する頂点 の確率密度の積であった。よって、サブパスの確率密度を乗算すればよいことになる。 – 例えば、eyeサブパスの頂点数3、lightサブパスの頂点数2なら、以下のような式の分母を計 算していることに相当する。 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓 𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝐺 𝑥1 ↔ 𝑥2 𝑓 𝑟 𝑥1 ← 𝑥2 ← 𝑦1 𝐺 𝑥2 ↔ 𝑦1 𝑓 𝑟 𝑥2 ← 𝑦1 ← 𝑦0 𝐺 𝑦1 ↔ 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 )𝑃𝐴 (𝑦0 )𝑃𝐴 (𝑦1 ) 光 源

Slide 147

Slide 147 text

頂点間の接続  続いて、モンテカルロ積分を解くためのスループットを求める。各サブパスのス ループットは既に求まっているため、サブパスの接続によるスループットを新た に計算すればよい。 – 例えば、eyeサブパスの頂点数3、lightサブパスの頂点数2なら、以下のような式の分子の黒い部分を 新たに計算する必要がある。今回はconnection_throughputが対応する。 – 基本的に、eyeサブパス側のBRDFと、lightサブパス側のBRDF、さらに端点間のG項で構成される。  ただし、lightサブパスの頂点数が1の時、非完全拡散光源の場合は相手の頂点の位 置次第でスループットが変化するため改めて光源からの放射輝度値を計算する。 – 今回は完全拡散光源なので単純にemissionの値を入れる。 𝑊′ 𝑥0 ← 𝑥1 𝐺 𝑥0 ↔ 𝑥1 𝑓 𝑟 𝑥0 ← 𝑥1 ← 𝑥2 𝐺 𝑥1 ↔ 𝑥2 𝑓 𝑟 𝑥1 ← 𝑥2 ← 𝑦1 𝐺 𝑥2 ↔ 𝑦1 𝑓 𝑟 𝑥2 ← 𝑦1 ← 𝑦0 𝐺 𝑦1 ↔ 𝑦0 𝐿𝑒 𝑦1 ← 𝑦0 𝑃𝐴 𝑥0 𝑃𝐴 𝑥1 𝑃𝐴 (𝑥2 )𝑃𝐴 (𝑦0 )𝑃𝐴 (𝑦1 ) 光 源

Slide 148

Slide 148 text

頂点間の接続  接続する際、サブパスの端点間に別の物体があれば当然そのパスのスルー プットは０になるので、端点間をレイトレーシングしておく。  ただし、lightサブパスを直接レンズにつなげるパターンというのがあるた め、実際に処理を打ち切るかどうかについてはもう少し後で判定する。

Slide 149

Slide 149 text

頂点間の接続  connection_throughput の計算を行っていく。まず、eyeサブパス側。  eyeサブパスの端点の材質によって処理が変わる。  まず、完全拡散面だった場合、 connection_throughputにBRDFを乗算する。 光 源

Slide 150

Slide 150 text

頂点間の接続  なお、端点が完全拡散面上だった場合かつ端点間に別の物体が存在すれば、 そこで処理を打ち切れる。 光 源

Slide 151

Slide 151 text

頂点間の接続  レンズ上だった場合、レンズへの入射方向によって最終的なイメージセン サ上の位置が決まるので、その計算を行う。また、 connection_throughput にW’を乗算する。 光 源

Slide 152

Slide 152 text

頂点間の接続  スペキュラ上において反射方向や屈折方向は決定的に決まる。この方向と、 lightサブパスの端点への方向が一致する確率は0なので、eyeサブパスの端 点がスペキュラ上に存在した場合、スループットは必ず0になる。（BRDFが 𝛿であり、0になるため）よって、そこで処理を打ち切る。  また、端点が光源上にいても打ち切る。（反射率0を仮定しているため） 光 源 反射方向 スペキュラ

Slide 153

Slide 153 text

頂点間の接続  次に、lightサブパス側について計算を行う。  eyeサブパス同様、lightサブパスの端点の材質によって処理が変わる。  まず、完全拡散面だった場合、 connection_throughputにBRDFを乗算する。 光 源

Slide 154

Slide 154 text

頂点間の接続  光源の反射率0を仮定しているため、ライトトレーシングの時、最初の頂点 以外は光源上に頂点生成されず、 lightサブパスの頂点数が1になる時以外で、 以下の赤枠部分に入ることは無い。 光 源

Slide 155

Slide 155 text

頂点間の接続  eyeサブパスの場合と同様に、サブパスの端点がスペキュラ上に存在した場 合、そこで処理を打ち切る。またレンズ上に存在した場合もスペキュラと みなすため、打ち切る。 光 源 スペキュラ 反射方向

Slide 156

Slide 156 text

モンテカルロ積分  最後に、端点間のG項を計算し、MISの重みを計算し、モンテカルロ積分を 行って終わり。  モンテカルロ積分は、今まで計算したスループットと端点間のスループッ トの積をパスのサンプリング確率密度で割ることによって行う。もちろん MISの重みも乗算する。  これにより、eyeサブパスの頂点数+lightサブパスの頂点数分の頂点数の積 分に対応したモンテカルロ積分が解かれたことになる。

Slide 157

Slide 157 text

MISの重み計算  あるeyeサブパスとlightサブパスから成るパスに対応するMISの重みを計算 することについて考える。 – lightサブパスの頂点数を𝑠、eyeサブパスの頂点数を𝑡とする。 – 𝑘 = 𝑠 + 𝑡 − 1とする。 – eyeサブパスを𝑥0 𝑥1 … 𝑥𝑡−1 、lightサブパスを𝑦𝑠−1 … 𝑦1 𝑦0 とする。 – 接続後のパスを𝑥 = 𝑥0 𝑥1 … 𝑥𝑡−1 𝑦𝑠−1 … 𝑦1 𝑦0 とする。  パワーヒューリスティックを使うと、重みは以下のようになる。 – ただし、𝑝𝑠,𝑡 𝑥 は頂点数𝑠 + 𝑡のモンテカルロ積分に対するサンプリング戦略の一つで、eye サブパスの頂点数𝑡、lightサブパスの頂点数𝑠の場合の𝑥 のサンプリング確率密度。 – 既に説明したとおり、eyeサブパスの頂点数𝑡 、lightサブパスの頂点数𝑠のモンテカルロ積分 は、頂点数𝑠 + 𝑡の積分に対するモンテカルロ積分の一つである（サンプリング戦略の一 つ）。この時、同じ頂点数𝑠 + 𝑡のモンテカルロ積分に対して、全部で𝑠 + 𝑡 + 1種類のサンプ リング戦略が考えられるので、以下のような式になる。 – 𝑠 + 𝑡 + 1種類のサンプリング戦略とは、すなわちeyeサブパス、lightサブパスの頂点数の組 み合わせの数である。 𝑤𝑠,𝑡 𝑥 = 𝑝𝑠,𝑡 𝑥 𝛽 𝑝𝑗, 𝑠+𝑡 −𝑗 𝑥 𝛽 𝑠+𝑡 𝑗=0

Slide 158

Slide 158 text

MISの重み計算の例  あるeyeサブパスとlightサブパスから成るパスに対応するMISの重みを計算 することについて考える。 – eyeサブパスの頂点数を𝟐、lightサブパスの頂点数を𝟑とする。  パワーヒューリスティックを使うと、重みは以下のようになる。 – 𝑥 を「仮に」eyeサブパスの頂点数5、lightサブパスの頂点数0の様なサンプリングによって得 られたとしたときの、確率密度が𝑝0,5 𝑥 。𝑝1,4 𝑥 や𝑝3,2 𝑥 等も同様。 𝑤3,2 𝑥 = 𝑝3,2 𝑥 𝛽 𝑝0,5 𝑥 𝛽 + 𝑝1,4 𝑥 𝛽 + 𝑝2,3 𝑥 𝛽 + 𝑝3,2 𝑥 𝛽 + 𝑝4,1 𝑥 𝛽 + 𝑝5,0 𝑥 𝛽

Slide 159

Slide 159 text

MISの重み計算の例  実際に得られたパスを、もし別のサンプリング戦略（頂点数の組み合わ せ）でサンプリングされたとしたら、どのような確率密度になるのかを計 算する必要があるということ。 光 源 実際のパス（eyeサブパス頂点数2、lightサブパス頂点数3） 光 源 同じパスを別の戦略（eyeサブパス頂点数4、lightサブパス頂点数1）でサンプリングした場合

Slide 160

Slide 160 text

calc_mis_weight()  以上の話に基づき、与えられたeyeサブパス、lightサブパスのMIS重みを計 算する。

Slide 161

Slide 161 text

calc_mis_weight()  まず、eyeサブパスとlightサブパスの始点である、光源上のサンプリング確 率密度とレンズ上のサンプリング確率密度を計算しておく。  さらに、後続の処理をやりやすくするためにeyeサブパスとlightサブパスを 一列に並べて置く。 – 先頭から順番に、lightサブパスの頂点を並べ、並べ終わったらeyeサブパスの頂点を逆順で 並べることで全体のパスとする。 – 𝑥 = 𝑥0 𝑥1 … 𝑥𝑡−1 𝑦𝑠−1 … 𝑦1 𝑦0 = 𝒛𝒔+𝒕−𝟏 … 𝒛𝟏 𝒛𝟎 – 𝑧𝑖 とvs[i]が対応している。

Slide 162

Slide 162 text

MISの重みの効率的計算  頂点数の合計が一定であることから、𝑝𝑠,𝑡 𝑥 を以下のように書き直す。 – 𝑖はあるサンプリング戦略における、lightサブパスの頂点数を示す。  𝑝𝑖 𝑥 を𝑖 = 0から𝑠 + 𝑡の範囲で計算すれば重みは計算できるが、単純にこれ を行うと効率が悪い。そこで、まず𝑝𝑖 𝑥 と𝑝𝑖+1 𝑥 の比、 𝑝𝑖+1 𝑥 𝑝𝑖 𝑥 を計算して、 その結果に基づいて𝑝𝑖 𝑥 を計算することにする。  ソースコード中ではこの比はpi1_piという配列が対応する。 𝑝𝑖 𝑥 = 𝑝𝑖, 𝑠+𝑡 −𝑖 𝑥

Slide 163

Slide 163 text

MISの重みの効率的計算  比の計算は以下のようにして行う。  𝑖 = 0の時、  ただし、𝑃𝐴 𝐸 𝑥 は、パストレーシングによって頂点𝑥がサンプリングされたと きの、面積測度に関する確率密度で、𝑃𝐴 𝐿 𝑥 は、ライトトレーシングによっ て頂点𝑥がサンプリングされたときの、面積測度に関する確率密度。  𝑃𝐴 𝐿 𝑦0 = 𝑃𝐴 𝐿(𝑧0 )は頂点が光源上からサンプリングされる確率密度で、これは 計算済み。  𝑃𝐴 𝐸 𝑥𝑘 = 𝑃𝐴 𝐸 𝑧0 を計算するわけだが、実際の分布は𝑃𝜎 𝐸(𝑧1 → 𝑧0 )に基づいて いるため、ひとつ前の頂点の情報も必要になる。さらに、𝑧1 がスペキュラ上 に存在した場合、さらにその一つ前の頂点𝑧2 の情報に必要になるので、都合 三つ分の頂点を使って𝑃𝐴 𝐸 𝑧0 を計算することになる。 – この計算はcalc_pdf_A()が行う。 – なお、頂点をlightサブパスの始点が最初に来るように一列に並べていたので、𝑥𝑘 は頂点配列上 では𝑧0 に対応し、𝑦0 も𝑧0 に対応している。 𝑝1 𝑥 𝑝0 𝑥 = 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−1 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−1 𝑃𝐴 𝐸(𝑥𝑘) = 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸(𝑥𝑘) = 𝑃𝐴 𝐿(𝑧0) 𝑃𝐴 𝐸(𝑧0)

Slide 164

Slide 164 text

MISの重みの効率的計算  なお、各頂点のサンプリング確率密度の計算に際してはロシアンルーレット による生成判定の確率も考慮する。 – ここで以下のように分母分子を打ち消しあうことが容易になるようにするために、ロシアン ルーレットの確率を頂点ごとに独立に計算していた。 𝑝1 𝑥 𝑝0 𝑥 = 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−1 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−1 𝑃𝐴 𝐸(𝑥𝑘) = 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸(𝑥𝑘) = 𝑃𝐴 𝐿(𝑧0) 𝑃𝐴 𝐸(𝑧0)

Slide 165

Slide 165 text

MISの重みの効率的計算  0 < 𝑖 < 𝑘の時、  𝑝𝑖 𝑥 はlightサブパス頂点数iのサンプリング戦略なので、同じパスについて、𝑦0 から𝑦𝑖−1 は𝑃𝐴 𝐿でサンプリングし、 𝑥𝑘−𝑖 から𝑥0 は𝑃𝐴 𝐸でサンプリングすることになる。  先ほどと同様に、対応する頂点のサンプリング確率密度を、その前の頂点の情報を使っ て計算していく。  一列に並べた頂点配列上で、𝑦𝑖 に対応するインデックスは𝑖で、その前の頂点は（ライト トレーシングなので）𝑖 − 1, 𝑖 − 2になる。つまり𝑧𝑖−2 , 𝑧𝑖−1 , 𝑧𝑖 。  𝑥𝑘−𝑖 に対応するインデックスも𝑖で、その前の頂点は（パストレーシングなので） 𝑖 + 1, 𝑖 + 2になる。つまり𝑧𝑖+2 , 𝑧𝑖+1 , 𝑧𝑖 。 𝑝𝑖+1 𝑥 𝑝𝑖 𝑥 = 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘− 𝑖+1 𝑃𝐴 𝐿 𝑦𝑖 𝑃𝐴 𝐿 𝑦𝑖−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−(𝑖+1) 𝑃𝐴 𝐸 𝑥𝑘−𝑖 𝑃𝐴 𝐿 𝑦𝑖−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) = 𝑃𝐴 𝐿(𝑦𝑖) 𝑃𝐴 𝐸(𝑥𝑘−𝑖) = 𝑃𝐴 𝐿(𝑧𝑖) 𝑃𝐴 𝐸(𝑧𝑖)

Slide 166

Slide 166 text

MISの重みの効率的計算  𝑖 = 𝑘の時、  𝑃𝐴 𝐸(𝑥0 ) は頂点がレンズ上からサンプリングされる確率密度で、これは計算済み。 𝑝𝑘+1 𝑥 𝑝𝑘 𝑥 = 𝑃𝐴 𝐿 𝑦𝑘 𝑃𝐴 𝐿 𝑦𝑘−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐿 𝑦𝑘−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) = 𝑃𝐴 𝐿(𝑦𝑘) 𝑃𝐴 𝐸(𝑥0) = 𝑃𝐴 𝐿(𝑧𝑘) 𝑃𝐴 𝐸(𝑧𝑘)

Slide 167

Slide 167 text

MISの重みの効率的計算  次に、比を使って𝑝𝑖 𝑥 を計算する。まず、初期値として𝑝𝑠 𝑥 を与える。こ れは、「実際に」サンプリングされたパスの確率密度なので、すでに求 まっている。  あとは、𝑝𝑖+1 𝑥 𝑝𝑖 𝑥 が全て分かっているため、 𝑝𝑠 𝑥 から初めてすべての𝑝𝑖 𝑥 を 計算することが出来る。

Slide 168

Slide 168 text

calc_pdf_A()  𝑃𝐴 𝐸 𝑥 と𝑃𝐴 𝐿 𝑦 を計算する関数。今回はパストレーシングもライトトレーシ ングも頂点のサンプリング自体はほとんど同じなので、共通の関数でよい。  基本的に、始点と次の頂点を与えられると次の頂点の面積測度に関するサ ンプリング確率密度を返す。スペキュラの考慮のため、始点のひとつ前の 頂点も与える。

Slide 169

Slide 169 text

calc_pdf_A()  始点の材質ごとに、次の頂点のサンプリング方法は変わる。始点が完全拡 散面や光源上に存在すれば、次の頂点はコサイン項によるインポータンス サンプリングによって行われる。  ここでの結果は、立体角測度に関するものなので、最後に面積測度に関す る確率密度に変換する。

Slide 170

Slide 170 text

calc_pdf_A()  始点がレンズ上に存在した場合、すでに見てきたようなやり方で次の頂点 のサンプリング確率密度を計算する。

Slide 171

Slide 171 text

calc_pdf_A()  始点がスペキュラ上に存在した場合、処理が分岐する。  まず、始点が完全鏡面上にいた場合は、次の頂点のサンプリング確率は𝛿𝜎 𝜔𝑖 − 𝜔𝑟 にな る。始点と次の頂点の組み合わせは、常に「実際に」サンプリングされたパスの一部で あるため、関数の値が0になることはない。 – 唯一、始点としてスペキュラ上のサブパスの端点が選ばれた場合のみ、次の頂点のサンプリング確率密度は0に なるはずだが、この場合はそもそもとっくに処理が打ち切られているはずなので考えなくて良い。  しかし、例によって数値的に処理しにくいので、𝛿𝜎 はあとでまとめて処理することにし て、ひとまず係数だけ考慮する。完全鏡面なら1になる。

Slide 172

Slide 172 text

calc_pdf_A()  始点がガラス上にいた場合、始点のひとつ前の頂点から始点へ向かうベクトル を考え、そのベクトルがガラスに進入しようとしているのか、出て行こうとし ているかを判定する。  これらの結果を使うことで、始点から次の頂点へのベクトルが反射によるもの なのか、屈折によるものなのかを判定することが出来る。そして、それぞれの イベントに応じてロシアンルーレットの確率を計算し、 𝛿𝜎 の係数とする。  また、全反射が起こる状況もある。この場合は𝛿𝜎 の係数は鏡面同様1になる。 from next prev_from

Slide 173

Slide 173 text

calc_pdf_A()  なお、ここの判定は、三頂点の位置関係から幾何的に計算しているが、パスト レーシング・ライトトレーシングの際に各頂点に屈折・反射・全反射といった 反射の種類に関するメタデータを付加しておけば、ここで改めて幾何的計算を 行う必要はない。  今回は、頂点に付加するメタデータが増えるのを嫌ってこのような実装にした が、これ以外の実装も当然考えられることに注意。 from next prev_from

Slide 174

Slide 174 text

スペキュラ処理  最後にスペキュラ処理をする。calc_pdf_A()の結果は、 𝛿𝜎 に関してはあくま で係数を保持しただけのものであった。ここで、一般に以下の式が成り 立っていた。 𝑝𝑖+1 𝑥 𝑝𝑖 𝑥 = 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘− 𝑖+1 𝑃𝐴 𝐿 𝑦𝑖 𝑃𝐴 𝐿 𝑦𝑖−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) 𝑃𝐴 𝐸 𝑥0 𝑃𝐴 𝐸 𝑥1 …𝑃𝐴 𝐸 𝑥𝑘−(𝑖+1) 𝑃𝐴 𝐸 𝑥𝑘−𝑖 𝑃𝐴 𝐿 𝑦𝑖−1 …𝑃𝐴 𝐿 𝑦1 𝑃𝐴 𝐿(𝑦0) = 𝑃𝐴 𝐿(𝑦𝑖) 𝑃𝐴 𝐸(𝑥𝑘−𝑖) = 𝑃𝐴 𝐿(𝑧𝑖) 𝑃𝐴 𝐸(𝑧𝑖)

Slide 175

Slide 175 text

スペキュラ処理  今、lightサブパスの端点𝑧𝑠−1 とeyeサブパスの端点𝑧𝑠 はスペキュラ上に存在 しない。このとき、各𝑝𝑖 について以下のような式が成り立つ。 – このとき、𝑝𝑠 𝑥 は実際にサンプリングされたパスの確率密度なのですでに求まっている。 – 𝐶(𝑧𝑠−1 → 𝑧𝑠 )は立体角測度を面積測度に変換するための例の係数（cos 𝜃 𝑟2 ） 𝑝𝑠+1 𝑥 = 𝑃𝐴 𝐿 𝑧𝑠 𝑃𝐴 𝐸 𝑧𝑠 𝑝𝑠 𝑥 = 𝑃 𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥 𝑝𝑠+2 𝑥 = 𝑃𝜎 𝐿 𝑧𝑠 → 𝑧𝑠+1 𝐶(𝑧𝑠 → 𝑧𝑠+1 ) 𝑃𝜎 𝐸 𝑧𝑠+2 → 𝑧𝑠+1 𝐶(𝑧𝑠+2 → 𝑧𝑠+1 ) 𝑃𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃 𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥 𝑝𝑠+3 𝑥 = 𝑃𝜎 𝐿 𝑧𝑠+1 → 𝑧𝑠+2 𝐶(𝑧𝑠+1 → 𝑧𝑠+2 ) 𝑃 𝜎 𝐸 𝑧𝑠+3 → 𝑧𝑠+2 𝐶(𝑧𝑠+3 → 𝑧𝑠+2 ) 𝑃𝜎 𝐿 𝑧𝑠 → 𝑧𝑠+1 𝐶(𝑧𝑠 → 𝑧𝑠+1 ) 𝑃𝜎 𝐸 𝑧𝑠+2 → 𝑧𝑠+1 𝐶(𝑧𝑠+2 → 𝑧𝑠+1 ) 𝑃 𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥

Slide 176

Slide 176 text

スペキュラ処理  仮に𝑧𝑠+1 がスペキュラ上に存在したとすると、𝑃𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 が𝛿𝜎 を含むこと になる。分子に𝛿𝜎 が存在しなければ、キャンセルできないので、分母に𝛿𝜎 が 残る𝑝𝑖 は0になる。  以下の場合、𝑝𝑠+1 と𝑝𝑠+2 が0になる。また、𝑝𝑠+3 は分母と分子に𝛿𝜎 が現れるの でキャンセルされ、0とは限らない。 𝑝𝑠+1 𝑥 = 𝑃𝐴 𝐿 𝑧𝑠 𝑃𝐴 𝐸 𝑧𝑠 𝑝𝑠 𝑥 = 𝑃 𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥 𝑝𝑠+2 𝑥 = 𝑃𝜎 𝐿 𝑧𝑠 → 𝑧𝑠+1 𝐶(𝑧𝑠 → 𝑧𝑠+1 ) 𝑃𝜎 𝐸 𝑧𝑠+2 → 𝑧𝑠+1 𝐶(𝑧𝑠+2 → 𝑧𝑠+1 ) 𝑃𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃 𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥 𝑝𝑠+3 𝑥 = 𝑃𝜎 𝐿 𝑧𝑠+1 → 𝑧𝑠+2 𝐶(𝑧𝑠+1 → 𝑧𝑠+2 ) 𝑃 𝜎 𝐸 𝑧𝑠+3 → 𝑧𝑠+2 𝐶(𝑧𝑠+3 → 𝑧𝑠+2 ) 𝑃𝜎 𝐿 𝑧𝑠 → 𝑧𝑠+1 𝐶(𝑧𝑠 → 𝑧𝑠+1 ) 𝑃𝜎 𝐸 𝑧𝑠+2 → 𝑧𝑠+1 𝐶(𝑧𝑠+2 → 𝑧𝑠+1 ) 𝑃 𝜎 𝐿 𝑧𝑠−1 → 𝑧𝑠 𝐶(𝑧𝑠−1 → 𝑧𝑠 ) 𝑃𝜎 𝐸 𝑧𝑠+1 → 𝑧𝑠 𝐶(𝑧𝑠+1 → 𝑧𝑠 ) 𝑝𝑠 𝑥

Slide 177

Slide 177 text

スペキュラ処理  今の話は𝑝𝑠+1 𝑥 , 𝑝𝑠+2 𝑥 , …についてだったが、𝑝𝑠−1 𝑥 , 𝑝𝑠−2 𝑥 , …について も同様。  まとめると、 頂点𝑧𝑖 がスペキュラ上に存在する場合、 𝑝𝑖 𝑥 , 𝑝𝑖+1 𝑥 は0にな るということ。これによって、 𝛿𝜎 を適切に処理したことになる。  また、直感的には、 𝑧𝑖 がスペキュラ上に存在すれば、その頂点を端点とす る接続によるパスのサンプリング確率は0になる、ということ。 𝑧𝑖 𝑧𝑖−1 𝑝𝑖 𝑥 𝑧𝑖 𝑧𝑖−1 𝑝𝑖+1 𝑥 いずれの場合も、 𝑧𝑖 がスペキュラ上に存在するなら、端点間を接続する 確率は0になる。

Slide 178

Slide 178 text

パワーヒューリスティック  最後に、以上の結果に基づいてパワーヒューリスティックによってMIS重み を計算する。

Slide 179

Slide 179 text

render.h  render_by_bidirectional_pathtracing()  双方向パストレーシングによるモンテカルロの結果を、具体的な画像に画 素値として与える。

Slide 180

Slide 180 text

render.h  得られたサンプルについて、サンプルが現在の画素 𝑥, 𝑦 から発射された eyeサブパスを含むものだった場合、 𝐼𝑥,𝑦 のモンテカルロ推定値は samples[i].valueそのものなので、そのまま足す。その後、下の画像出力時に 発射された回数の総計(iteration_per_thread * num_threads)で割る。  得られたサンプルについて、現在の画素から発射されたeyeサブパスを含む ものではなかった場合（lightサブパスが別の画素 𝑥′, 𝑦′ に到達した場合）は 𝐼𝑥′,𝑦′ のモンテカルロ推定値を新しく得たわけだが、この場合、画像全体に 対して光源側からサンプルを生成し、たまたま 𝑥′, 𝑦′ 'にヒットしたと考え るため、このようなサンプルについては最終的に光源から発射した回数の 総計(width * height * iteration_per_thread * num_threads)で割って、画素への 寄与とする必要がある。  iteration_per_thread * num_threadsの分は上と共通なので、width * heightで 割ってからimage_bufferに足すことで、最終的な画像出力時に帳尻があり、 正確な結果になる。

Slide 181

Slide 181 text

7. レンダリング結果

Slide 182

Slide 182 text

比較 パストレーシング 59秒 ライトトレーシング 62秒 双方向パストレーシング 60秒 コーネルボックス

Slide 183

Slide 183 text

比較 小光源 パストレーシング 61秒 ライトトレーシング 61秒 双方向パストレーシング 63秒

Slide 184

Slide 184 text

比較 ガラス球 パストレーシング 61秒 ライトトレーシング 66秒 双方向パストレーシング 67秒

Slide 185

Slide 185 text

収束 1000 10 250 50 コーネルボックス 双方向パストレーシングの実行回数/pixel

Slide 186

Slide 186 text

謝辞  アドバイスや感想を下さった以下の方々に感謝いたします。 – Yuki OZAWA（@qatnonoil） – お餅（@omochi64）

Slide 187

Slide 187 text

参考文献  Eric Veach “Robust Monte Carlo Methods for Light Transport Simulation “  Dietger van Antwerpen “A Survey of Importance Sampling Applications in Unbiased Physically Based Rendering”  Philip Dutre, Philippe Bekaert, Kavita Bala “Advanced Global Illumination, Second Edition”  Matt Pharr, Greg Humphreys “Physically Based Rendering, Second Edition: From Theory To Implementation”  Kevin Suffern “Ray Tracing from the Ground Up”  Computer Graphics “http://www35.atpages.jp/shocker/memoRANDOM/CG/CG.php”