本記事ではパストレーシングの基本的なノイズ低減手法について説明します．なお，前の記事を理解している前提の説明になります．
前の記事で説明したパストレーシングでは，入射方向のサンプリングとヒット位置の計算を繰り返して，運よくライトにヒットすることで，ようやくライトの影響が計算できました．しかしながらこの方法では，下の図のようにライトが非常に小さかったりライトが遠くにあると，ライトにヒットする確率が非常に小さくなり，ライトからカメラに至る光の経路をほとんどサンプリングされないという問題があります．
この問題を解決するための手法が本記事で説明するNEE（Next Event Estimation）です．この方法を用いることで，ライトが小さかったり，ライトが遠くにあっても安定してライトの影響を計算できます．
しかしながらこの手法では，逆にライトが大きかったり，ライトが近くにあるとあまりいい効果が得られず，材質によっては偶然のヒットに頼るよりも圧倒的に悪い結果になってしまうこともあります．
そこで使用する手法が多重重点的サンプリング（MIS：Multiple Importance Sampling）です．この手法を用いることでNEEと偶然のヒットに頼る方法をうまく組み合わせることができ，安定した推定が可能になります．
本記事では，NEEを説明したのちにMISの説明をし，最後にパストレーシングの発展形である双方向パストレーシングの導入も説明します．
NEE: Next Event Estimation
NEEでは，サンプリングした入射方向が偶然ライトにヒットすることに頼るのではなく，ライト上の位置をサンプリングし輝度計算点と明示的に接続してしまいます．
ライト上の位置を直接サンプリングしてしまうので，ライトが小さかったり，ライトが遠くにあっても安定して影響を計算することができます．
アイディアは以上なのですが，ここでサンプリングするのは発光物体上の位置であり，方向ではないので，方向に関する積分であるレンダリング方程式をそのままでは使用できません．
そのため，レンダリング方程式を位置に関する積分に書き換える必要があります．立体角は単位球上の面積だったので，下図に示すように，微小立体角と微小面積には以下の関係があります．
\begin{aligned}
d\omega &= \frac{\cos\theta'}{||x-y||^2}dA
\end{aligned}
同じ面積の面であっても，遠くにあるほど，また，その面がこちらを向く方向から傾いていくほど見かけの面積，つまりは立体角が小さくなるということを表しています．この関係式を使うことで，レンダリング方程式の反射輝度は以下のように変形できます．
\begin{aligned}
L_r(x,\omega_o) &= \int_\Omega L_i(x,\omega_i)f(x,\omega_o,\omega_i)\cos\theta d\omega_i \\
&= \int_\mathcal{M} L_o(y,\omega_{yx})f(x,\omega_o,\omega_{xy})\frac{\cos\theta\cos\theta'}{||x-y||^2}V(x,y) dA \\
&= \int_\mathcal{M} L_o(y,\omega_{yx})f(x,\omega_o,\omega_{xy})G(x,y)V(x,y) dA \\
\end{aligned}
ここで，$\mathcal{M}$はシーン上の物体表面，$L_i$は入射輝度，$L_o$は出射輝度，$f$はBRDF，$G$は幾何項，$V$は可視関数，$\omega_{yx}$は位置$y$から位置$x$への方向です．
可視関数$V(x,y)$は2点間の間に遮蔽物があれば0，それ以外なら1となる関数です．位置に関する積分では輝度計算点$x$とシーン上の点$y$の間に遮蔽物がある場合があるため，そのような位置$y$からの影響を消すために可視関数が必要になります．
上の図では黒太線の面は位置$x$から見えていないので可視関数$V$は0となります．$V=1$の微小面を位置$x$を中心とする単位球に投影していけば，ちょうど上半球が埋まり，位置に関する積分も方向に関する積分も同じものを表していそうです．
NEEは光源上の位置との接続だったので，反射輝度は以下の式となります．
\begin{aligned}
L_r(x,\omega_o,\omega_i) &= \int_\Omega \left\{ L_e(r(x,\omega_i),-\omega_i)+ L_r(r(x,\omega_i),-\omega_i) \right\} f(x,\omega_o,\omega_i)\cos\theta d\omega_i \\
&= \int_\mathcal{M_e} L_e(y,\omega_{yx})f(x,\omega_o,\omega_{xy})G(x,y)V(x,y) dA \\
&+ \int_\Omega L_r(r(x,\omega_i),-\omega_i) f(x,\omega_o,\omega_i)\cos\theta d\omega_i \\
\end{aligned}
ここで，$L_e$は発光輝度，$r(x,\omega_i)$は位置$x$から方向$\omega_i$にレイを飛ばしたときのヒット位置を返す関数です．
この式で重要なのは，方向に関する積分ではヒット位置で発光輝度の計算をしないということです．これは発光輝度の成分はNEEで計算しているためです．どちらでも計算してしまうと，発光輝度の影響が2倍になってしまいます．
NEEで解く積分値を$L_d$として，モンテカルロ法を用いたときの推定式は以下のようになります．なお，$L_d$の$d$は直接照明（Direct Illumination）の$d$です．
\begin{aligned}
L_d\approx\hat{L}_d= \frac{1}{N_d}\sum_{s=1}^{N_d}\frac{L_e(y_s,\omega_{yx})f(x,\omega_o,\omega_{xy})G(x,y)V(x,y)}{p(y_s|x,\omega_o)}
\end{aligned}
ここで，$N_d$はサンプル数，$y_s$は$s$回目の試行における位置のサンプルで，確率密度関数$p$に従います．
確率密度関数は被積分関数に比例したものを使用するのが理想でした．特定の条件下ではそれが可能ですが，一般的には，面積について一様な確率でサンプリングしたり，ライトのエネルギーに比例した確率密度関数を使用します．
ライトが三角形の集合で表現されている場合は，まず三角形をサンプリングし，その三角形上を一様にサンプリングすることで位置を決定します．確率密度は（三角形を選択する確率）×（三角形上の点を選択する確率密度）です．
問題点
上で説明したように，NEEを使えば小さいライトがあっても，遠くにライトがあっても安定した推定ができそうです．しかしながら，NEEを使うことによって結果が悪くなるケースがあります．それは，ライトが近くにある場合や，鋭い反射をする材質の場合です。これらはNEEを使わなかった時には苦手ではなかったケースです。なぜ，NEEを使うことによってこれらのケースで結果が悪くなってしまうのでしょうか？
まずライトが近くにある場合です．下の図のように，輝度計算点$x$の近くに発光する面があると，その面が占める立体角が支配的となり，反射輝度に大きく寄与するようになります．しかしながらNEEでは，輝度計算点$x$の近くの位置を重点的にサンプリングするということはしないため，確率密度関数が被積分関数に比例するほどよいという理想から大きく外れてしまい推定は安定しません．
ほとんどは値が小さく推定され，たまに輝度計算点$x$に近い位置$y$がサンプリングされたときに，幾何項に含まれる$1/||x-y||^2$によって異常に明るく推定されます．入射方向サンプリングで推定していたときは，立体角が大きいライトは自然とヒットする確率が高くなるので問題がありませんでした．
次に鋭い反射をする材質の場合です．下の図で緑の楕円のようなものはBRDFの入射方向による値の大きさを表していて，鏡のような鋭い反射をする材質では，正反射方向付近の狭い範囲の光しか出射方向$\omega_o$に影響を与えません．しかしながら先ほどと同様，NEEではBRDFを考慮した位置のサンプリングを行わないため，確率密度関数が被積分関数に比例するほどよいという理想から大きく外れてしまします．
ほとんどは値が小さく推定され，たまにBRDFの値が大きい方向に位置がサンプリングされると，BRDFの値によって異常に大きく推定されます（鋭い反射をする材質のBRDFは非常に大きい値をとります）．入射方向サンプリングで推定していたときは，BRDFの値が大きい方向を重点的にサンプリングしていたため，BRDFが大きいということは確率密度関数値も大きくなり，分母分子で相殺されて適度な値に落ち着くため問題がありませんでした．
さて，一見よさそうだったNEEですがいろいろと問題がることが分かりました．そして，NEEが苦手とするケースは入射方向サンプリングが得意で，NEEが得意とするケースは入射方向サンプリングでは苦手でした．
そうすると考えることと言えば，両方を使っていい感じにしたい！ということでしょう．しかしながら先に説明したように，単に両方を使うだけでは，発光輝度の影響を2重に計算してしまいます．
では，両方を使ったうえでそれぞれ0.5倍にするのはどうでしょうか？たしかに2重に計算していた問題は解決されますが，両方の利点を活かしているとはいいがたいです．各手法は苦手なケースでは推定値が真値の数万倍になることもあり，それが0.5倍になったところで高精度な推定は期待できません．この問題を解決するのが次に説明する多重重点的サンプリングです．
多重重点的サンプリング（MIS:Multiple Importance Sampling）
MISは複数のサンプリング方法を使用して積分値や総和を推定する手法です．なおMISにおいては，サンプリング方法のことをサンプリング戦略と呼ぶのが一般的なので，ここでもサンプリング戦略と呼びます．
MISでは，各戦略は被積分関数の積分の推定を行うのではなく，被積分関数と重み関数の積の積分を推定します．最終的な積分値の推定値は各戦略の推定値の和となります．
\begin{aligned}
I&= \sum_{i=1}^{N}\int w_i(x)f(x)dx \\
I&\approx\hat{I}= \sum_{i=1}^{N}\frac{1}{n_i} \sum_{j=1}^{n_i}w_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})}
\end{aligned}
ここで，$N$は戦略数，$n_i$は$i$番目の戦略のサンプル数，$w_i$は$i$番目の戦略に対する重み関数，$X_{i,j}$は$i$番目の戦略の$j$番目のサンプルで，確率密度関数$p_i$に従います．
以下に示すように，重み関数の総和が$1$であれば，このように分解しても結果は元の積分値と一致します．
\begin{aligned}
I &= \int \left\{\sum_{i=1}^{N} w_i(x) \right\} f(x)dx \\
&= \int f(x)dx
\end{aligned}
イメージ図は以下の通りです．分かりやすさのため$f(x)=$定数としています．
通常のモンテカルロ法では，被積分関数$f(x)>0$ならば確率密度関数$p(x)>0$である必要がありました．$p(x)=0$にもかかわらず$f(x)>0$の領域があると，そこの影響を見逃してしまうからです．
一方MISでは複数の戦略があるため，いずれかの戦略で$f(x)>0$ならば$p_i(x)>0$を満たせば十分です．なお，$p_i(x)=0$ならば$w_i(x)=0$である必要があります．さもないと，$w_i(x)f(x)$の影響を見逃してしまうからです．
重み関数の導出
さて，MISでは各戦略に一様な重みを割り当てるのではなく，$x$ごとに重みを割り当てます．これにより，各戦略をより柔軟に組み合わせることが可能になりましたが，どのような重み関数を用いるのが最適なのでしょうか？
推定の目標は真値を高精度に推定することです．つまり，推定値の分散を小さくするのが目標となります．どのような重み関数が分散を小さくするのでしょうか？MISの推定値の分散は以下の式で表されます．
\begin{aligned}
V \left[ \hat{I} \right] &= V \left[ \sum_{i=1}^{N}\frac{1}{n_i} \sum_{j=1}^{n_i}w_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})} \right] \\
&= \sum_{i=1}^{N}\frac{1}{n_i} V \left[ w_i(X_{i})\frac{f(X_{i})}{p_i(X_{i})} \right] \\
&= \sum_{i=1}^{N}\frac{1}{n_i} \left\{ E \left[ w_i(X_{i})^2\frac{f(X_{i})^2}{p_i(X_{i})^2} \right] - E \left[ w_i(X_{i})\frac{f(X_{i})}{p_i(X_{i})} \right]^2 \right\} \\
\end{aligned}
この分散を最小化する重み関数を導出したいのですが，このままでは困難であるため，分散の上限値を最小化することにします．期待値の2乗の項は常に正であるため，分散の上限値は以下のようになります．
\begin{aligned}
V \left[ \hat{I} \right] &\le \sum_{i=1}^{N}\frac{1}{n_i} E \left[ w_i(X_{i})^2\frac{f(X_{i})^2}{p_i(X_{i})^2} \right] \\
&= \sum_{i=1}^{N}\frac{1}{n_i} \int_D w_i(x)^2\frac{f(x)^2}{p_i(x)^2}p(x)dx \\
&= \sum_{i=1}^{N}\frac{1}{n_i} \int_D w_i(x)^2\frac{f(x)^2}{p_i(x)}dx \\
\end{aligned}
分散の上限値を最小化するにあたって，ある点$x$での重み関数の値は，他の点$x'$に影響を与えないため，最終的に最小化するべき式は以下のようになります．
\begin{aligned}
\sum_{i=1}^{N} w_i(x)^2\frac{f(x)^2}{n_i p_i(x)} \\
\end{aligned}
この式を，重み関数の総和の制約（$\sum_{i=1}^Nw_i(x)=1$）の下で最小化すれば，最適ではないですが，準最適な重み関数が得られます．なお，他にも$p_i(x)=0$ならば$w_i(x)=0$などの制約がありましたが，総和の制約を満たすように導出した重み関数は他の制約も勝手に満たしてくれています．
実際にこの式を最小化する重み関数を導出したいと思います．これにはラグランジュ未定乗数法という方法を使います．自分はなぜこれで解けるのかの分かりやすい説明ができないので，手順だけ書きます．以下のように最小化対象の式に（制約関数×謎の変数$\lambda$）を加算した関数$O$を作ります．
\begin{aligned}
O(w_i)&=\sum_{i=1}^{N} w_i^2\frac{f^2}{n_i p_i} + \lambda \left( 1-\sum_{i=1}^N w_i \right)
\end{aligned}
そして，関数$O$を$w_i$と$\lambda$でそれぞれ偏微分して0と置いて，連立方程式を解くと目的の$w_i$が導出できます．まず，$w_i$についてです．
\begin{aligned}
\frac{\partial O}{\partial w_i} = \frac{2w_if^2}{n_ip_i}-\lambda =0\\
w_i = \frac{\lambda n_ip_i}{2f^2}
\end{aligned}
つぎに$\lambda$についてです．
\begin{aligned}
\frac{\partial O}{\partial \lambda} = 1-\sum_{i=1}^N w_i = 0 \\
\sum_{i=1}^N \frac{\lambda n_ip_i}{2f^2} = 1 \\
\lambda=\frac{2f^2}{ \sum_{i=1}^N n_ip_i} \\
\end{aligned}
さきほどの$w_i$の式に$\lambda$を代入して
\begin{aligned}
w_i &= \frac{n_ip_i}{2f^2} \frac{2f^2}{ \sum_{i=1}^N n_ip_i} \\
&= \frac{n_ip_i}{ \sum_{i=1}^N n_ip_i} \\
\end{aligned}
これが分散の上限値を最小化する重み関数となります．ちなみに，このようにして導出された重み関数をBalance Heuristicといいます．$x$をサンプリングしやすい戦略ほど，また，サンプリング回数が多い戦略ほど重みが大きくなります．
以下の図は，この重み関数を用いたときの各戦略が解く関数のイメージです．
各戦略は得意不得意があり，単独では不得意な領域のせいで推定精度が悪くなりますが，重み関数のおかげで，得意な領域に集中することができています．実際にサンプリングを行った場合は以下のようになり，安定した推定ができるのが分かります．
反射輝度の推定にMISを用いた場合の式は以下のとおりです。
\begin{aligned}
L_r(x,\omega_o,\omega_i) &= \int_\mathcal{M_e} w_{NEE}(y) L_e(y,\omega_{yx})f(x,\omega_o,\omega_{xy})G(x,y)V(x,y) dA \\
&+ \int_\Omega w_{BRDF}(\omega_i) L_e(r(x,\omega_i),-\omega_i) f(x,\omega_o,\omega_i)\cos\theta d\omega_i \\
&+ \int_\Omega L_r(r(x,\omega_i),-\omega_i) f(x,\omega_o,\omega_i)\cos\theta d\omega_i \\
\end{aligned}
ここで，$w_{NEE}$と$w_{BRDF}$はそれぞれNEEとBRDFサンプリングにおけるMIS重み関数です．実際にモンテカルロ法で推定する際は，2つ目と3つ目の積分の推定に使う入射方向サンプルは共有しても問題ないため，NEEだけを使ったときと計算コストは大きく変わらないです．
重み関数を計算するうえで注意が必要なのは，重み関数で使う確率密度関数は位置なら位置、方向なら方向に関する確率密度で揃えておく必要があるということです．位置に揃えるなら以下の関係式を用いて方向に関する確率密度を変換します．
\begin{aligned}
p(y)dA&=p(\omega)d\omega \\
p(y)&=p(\omega) \left| \frac{d\omega}{dA} \right| \\
&=p(\omega) \frac{\cos\theta'}{||x-y||^2} \\
\end{aligned}
双方向パストレーシングの導入
双方向パストレーシングは，その名の通りカメラ側からだけ経路をサンプリングするのではなく，ライト側からも経路をサンプリングする手法です．
実はNEEもライト側から経路の起点となる最初の位置をサンプリングしているので，双方向パストレーシングの一種と言えなくもないです．NEEによって推定の安定性が大きく向上しますが，それでもまだまだ十分ではありません．例えば以下の図のようにライトが狭い範囲を照らしている場合を考えます．
位置$x_2$がサンプリングされたとき，確かにNEEによって$x_2$での反射輝度を高精度に推定できます．
しかしながら，そもそも位置$x_1$においてライトの影響を受ける位置をサンプリングするには赤破線の間の狭い範囲の方向をサンプリングする必要があり，その確率は非常に小さくなります．つまり，肝心の$x_1$での反射輝度の推定精度は低いままです．NEEはライトからの直接照明には効果がありますが，光が反射して届く間接照明の改善はできないのです．
では，ライトから1回だけレイを飛ばして，NEEの要領でそのヒット位置$y_1$と位置$x_1$を接続してみるのはどうでしょうか？
位置$x_1$においてライトが照らす狭い範囲の方向をサンプリングする必要はなくなり，位置$x_1$に入射する間接照明の反射輝度を高精度に推定できそうです．
ただし注意なのは，この方法の推定精度が常によりよいとは限らず，今まで通りNEE＋入射方向サンプリングの方がいいこともあります．双方向パストレーシングでは，どれがいい方法かをがんばって何とか判定するなんてことはせず，MISで勝手にいい感じに重みを付けて対応します．
双方向パストレーシングでは，上のパターンだけでなく，ライトからより長い経路をサンプリングしたパターンも用います．通常のパストレーシングより使用しているサンプリング戦略の数が圧倒的に多いので，明らかに双方向パストレーシングの方が優秀そうです．説明は別の記事にします．