23. モンテカルロレンダリング

第12章ではパストレーシング法に出会いました：目からランダムな経路を飛ばして平均すれば、レンダリング方程式が解ける。
但し、ノイズは$1/\sqrt{N}$でしか減らない。
其の章は其れが上手く行く事を見せ、詳しい話は第2部で、と約束しました。
此の章が其の約束です。

モンテカルロ積分を基礎から組み立て（確率変数、期待値、分散、推定量）、次に全ての本格的なレンダラーが少ないサンプルで綺麗な画像を得る為に使う道具を見ます：逆関数法によるサンプリング、重点的サンプリング（コサイン、GGX、光源）、層化サンプリングと準モンテカルロ法、そして多重重点的サンプリングです。
其の後、パストレーシングを経路の積分として見直し、スループット、ロシアンルーレット、鏡、そして小さなC++のパストレーサーを扱います。
更に先へ進む手法（双方向、メトロポリス、プログレッシブフォトンマッピング）、バイアスと一致性、最後にデノイズと適応的サンプリングです。

第2部は検定の出題範囲を超えますが、此の章は其の項目の1つ、大域照明計算を深める物です：パストレーシング法、モンテカルロ積分、重点的サンプリング、ロシアンルーレットは全て第12章の検定のポイントに出てきました。
最初の幾つかの節で、其れらが何故検定で問われる通りに振る舞うのかを説明します。
前提は第12章のレンダリング方程式、BRDF、光路の表記で、23.5節では第22章のGGX分布を使います。
数学は学校の確率と第12章の積分です。
積分は必ず言葉でも読み上げます。

23.1 何故レンダリングに乱数が必要なのか

第12章（12.7節）のレンダリング方程式を思い出して下さい：

1つのピクセルは、此の式を1回評価すれば済む訳ではありません。
ピクセルの値は、ピクセルの面積（アンチエイリアシング）、レンズ（被写界深度）、シャッターが開いている時間（モーションブラー）についての積分で、更に跳ね返る度に入射方向の半球についての積分、ソフトシャドウの為にはランプの面積についての積分も加わります。
其の1つ1つが次元を増やします：ピクセルの位置で2、レンズで2、時間で1、そして跳ね返り1回毎に新しい方向で2、ランプ上の点で2。
3回跳ね返ると

各次元に10点の格子を取ると${10}^{17}$回の評価が必要で、1秒に10億本のレイを追跡するGPUでも約3年、其れも1ピクセル分です。
更に悪い事に、台形公式の誤差は1次元なら$N^{-2}$で減りますが、$d$次元では$N$個の点が$d$本の軸に分けられる為、$N^{-2/d}$でしか減りません。
$d=17$なら$N^{-0.12}$で、仕事を2倍にしても結果は9%も良くなりません。
此れが次元の呪いです。

モンテカルロ積分は此の呪いを解きます。
其の誤差は次元に関係なく$N^{-1/2}$で減り、サンプル数は幾つでも良く、高次の格子公式が強みを失う被積分関数の不連続（影の境界、物体の輪郭）も気にしません。
4次元を超えれば台形公式に勝ち、レンダリングは4次元より遥かに上にあります。
名前は1940年代のロスアラモスに由来します。
スタニスワフ・ウラム、ジョン・フォン・ノイマン、ニコラス・メトロポリスが中性子の輸送をランダムなサンプリングで計算し（Metropolis and Ulam, "The Monte Carlo Method", 1949）、其れは光の輸送の近い親戚です。

23.2 確率の基本

確率から必要なのは、5つの考え方だけです。

確率変数と期待値

確率変数とは、ランダムな試行から出てくる数です：サイコロの目$X$や、パストレーサーが選ぶ方向です。
其の期待値$E[X]$は、無限回繰り返した時の平均です。
公平なサイコロなら

分散と標準偏差

分散は、値が平均の周りにどれ位散らばるかを測ります：平均からの距離の2乗の期待値です。

サイコロでは$E[X^2]=(1+4+9+16+25+36)/6=91/6=15.1667$なので、$\mathrm{Var}[X]=15.1667-{3.5}^2=2.9167$です。
其の平方根が標準偏差で、$\sigma =1.708$、$X$と同じ単位を持ちます。
画像では、ピクセルの推定値の標準偏差こそが、目に見えるノイズです。
（第12章の光の「分散」と同じ言葉ですが、別の意味です。）

連続な確率変数

方向や位置は、或る範囲の中のどんな値にもなり得ます。
其の確率は確率密度関数、pdf $p(x)$で表します：$X$が$x$の周りの小さな区間$\mathrm{d}x$に入る確率が$p(x)\mathrm{d}x$です。
pdfは負にならず、積分すると1になり、1より大きくなる事もあります（幅0.1の区間に押し込められたpdfの高さは10です）。
左端から$x$迄の積分が累積分布関数、CDF $P(x)$で、$X\le x$となる確率です。

和は積分になります：

一番簡単なのは$[0,1)$上の一様分布$p(x)=1$で、コンピュータの乱数生成器がくれるのは此れです。
平均は$1/2$、分散は${\int }_0^1{x}^2\mathrm{d}x-(1/2{)}^2=1/3-1/4=1/12$です。
此の章では$U$（具体的な値なら$u$）は常に此の一様乱数を表します。

3つの規則

誤差の解析に必要なのは、次の3つの規則だけです：

• 期待値は線型：$E[aX+bY]=aE[X]+bE[Y]$。常に成り立つ。
• 定数倍すると分散は2乗倍：$\mathrm{Var}[aX]=a^2\mathrm{Var}[X]$。
• 独立な確率変数なら分散は足し算：$\mathrm{Var}[X+Y]=\mathrm{Var}[X]+\mathrm{Var}[Y]$。

従って、分散${\sigma }^2$の確率変数の独立なコピーを$N$個平均すると

平均は其の儘で、散らばりは$\sqrt{N}$で縮みます。
此れが第12章の$1/\sqrt{N}$の法則で、2番目の規則の2乗から来ています。

23.3 モンテカルロ推定量

$I=\int f(x)\mathrm{d}x$を推定するには、pdf $p$で独立なサンプル$X_k$を$N$個取り、次を計算します：

$F_N$は推定量です：其の値が$I$の推測になる確率変数です。
$p$で割るのが肝心な所です：サンプルが密な所（$p$が大きい）では1つ1つの重みが小さく、疎な所では大きくなります。

何故平均すると正しいのか

1つの項の期待値は

$p$が打ち消し合い、線型性から$N$個の項の平均も同じ期待値を持ちます。
期待値が丁度真の値になる推定量を不偏と言います。
打ち消しには条件が1つ隠れています：$f(x)\ne 0$となる所では必ず$p(x)>0$。
一度もサンプルされない領域は、結果から黙って抜け落ちます。

1回の推定はどれだけずれるか

23.2節の規則から、$F_N$の分散は${\sigma }^2/N$です。
${\sigma }^2$は1つの項の分散です：

此の章の残りでノイズを減らす為にする事は全て、${\sigma }^2$を小さくする（より良い$p$、23.5節、23.7節、23.8節）、サンプルを上手く独立でなくする（23.6節）、又は後からノイズを取り除く（23.12節）の何れかです。

計算例：正弦の山の面積

積分

は答えが分かっているので、推定量の働きを観察出来ます。
$X=\pi U$と$[0,\pi ]$で一様に選ぶと$p(x)=1/\pi$で、各項は$\pi \sin X$です。
4つの乱数$u=0.62,0.35,0.48,0.81$では：

33%大き過ぎます：4つの内3つのサンプルが山の頂上付近に落ちたのです。
期待される散らばりは

（${\int }_0^{\pi }{\sin }^{2}x\mathrm{d}x=\pi /2$を使いました）。
$N=4$なら$F_4$の標準偏差は$0.9669/2=0.483$なので、0.655の誤差は運が悪いものの珍しくはありません：標準偏差の約1.4倍です。

$N$が大きければ推定値は釣鐘型の分布に従い（中心極限定理）、約95%が真の値の$\pm 1.96\sigma /\sqrt{N}$以内に入ります：正弦で$N=100$なら$\pm 0.19$です。
レンダラーは$\sigma$を知りませんが、サンプル自身から推定出来ます。
適応的サンプリング（23.12節）が使うのは其れです。

検定のポイント

此の節はシラバスの大域照明計算の項目（第12章12.10節）に繋がります：

• モンテカルロ推定量は$f/p$を平均する。$f\ne 0$の所で$p>0$なら不偏。
• 標準偏差は$\sigma /\sqrt{N}$。サンプルを4倍にするとノイズは半分、100倍で正しい桁が1つ増える。
• $N^{-1/2}$という速さは次元に依存しない。其の為、大域照明の高次元の積分にモンテカルロ法が使われる。

23.4 乱数からサンプルを作る

コンピュータがくれるのは$[0,1)$の一様な数$u$です。
他の分布は全て其処から作らなければなりません。

逆関数法

主力は逆関数法です：CDF $P$を求め、次の様にします。

何故上手く行くのでしょうか？
$X\le x$となるのは丁度$U\le P(x)$の時で、一様な$U$なら其の確率は$P(x)$です。
従って$X$のCDFは$P$、つまりpdfは$p$です。
手順はいつも同じです：pdfを積分して$P$を求め、$P(x)=u$と置き、$x$について解きます。

例1：霧の中の距離。
密度$\sigma$の媒質の中で、光子が相互作用する迄に飛ぶ距離$t$のpdfは$p(t)=\sigma e^{-\sigma t}$です（第12章12.6節のランベルト・ベールの法則を確率として読んだ物）。
$P(t)=1-e^{-\sigma t}$なので、$1-e^{-\sigma t}=u$を解くと

1メートルあたり$\sigma =0.5$、$u=0.5$なら$t=\ln 2/0.5=1.386$ m、自由行程の中央値です。
第28章のボリュームのパストレーサーは、正に此の1行を使います。

例2：三角形。
後で$[0,\pi ]$上の三角形のpdfが必要になります。
$\pi /2$の頂点迄直線的に上がり、又下がる物です：

（三角形の底辺は$\pi$、高さは$2/\pi$なので面積は1です。）
左半分では$P(x)=2x^2/{\pi }^2$で、頂点で$1/2$になります。
$2x^2/{\pi }^2=u$を解くと、では$x=\frac{\pi }{2}\sqrt{2u}$、対称性から$u\ge 1/2$では$x=\pi -\frac{\pi }{2}\sqrt{2(1-u)}$です。

$P$の逆関数が求まらない時は、棄却法が代わりになります：形を囲む簡単な領域で候補を引き、外れた物を捨てるのです。
単位円板の中の点なら、正方形$[-1,1{]}^2$で$(x,y)$を引き、$x^2+y^2\le 1$なら採用します。
円板は正方形の$\pi /4=78.5\%$を占めるので、1点あたり平均1.27個の候補が要ります。
ループの回数が一定しない事が、GPUが逆関数法を好む理由です。

半球上の方向

法線の周りの半球上の方向は、極角$\theta$（法線から測る）と方位角$\varphi$（法線の周り）で表します。
方向の小さなパッチの立体角は$\mathrm{d}\omega =\sin \theta \mathrm{d}\theta \mathrm{d}\varphi$です：地平線近くの輪は、極の近くの輪より長いのです。
立体角についてのpdf $p(\omega )$は、$(\theta ,\varphi )$についてはpdf $p(\omega )\sin \theta$になり、$\theta$と$\varphi$を別々に逆関数法でサンプリング出来ます。

半球上の一様サンプリング：$p(\omega )=1/(2\pi )$。
$\theta$のpdfは$\sin \theta$、CDFは$1-\cos \theta$なので$\cos \theta =1-u_1$（$1-u_1$も一様なので、単に$\cos \theta =u_1$でも良い）、$\varphi =2\pi u_2$です。
つまり$z=\cos \theta$は$[0,1]$で一様で、第12章が使った事実です。

コサイン重み付き半球サンプリング：$p(\omega )=\cos \theta /\pi$（${\int }_{\mathrm{\Omega }}\cos \theta \mathrm{d}\omega =\pi$なので、$\pi$で正規化します）。
$\theta$のpdfは$\frac{\cos \theta }{\pi }\sin \theta \cdot 2\pi =2\sin \theta \cos \theta$、CDFは${\sin }^2\theta$なので

局所座標（法線がz軸）では$(\sqrt{u_1}\cos \varphi ,\sqrt{u_1}\sin \varphi ,\sqrt{1-u_1})$です：単位円板上の一様な点を、真上に持ち上げて半球に乗せた物です（マレーの方法）。

計算例：$u_1=0.25$、$u_2=0.5$なら$\sin \theta =0.5$、$\cos \theta =\sqrt{0.75}=0.8660$（$\theta =30\mathrm{°}$）、$\varphi =\pi$で、方向は$(-0.5,0,0.8660)$、pdfは1ステラジアンあたり$0.8660/\pi =0.2757$です。

シーンの中では、法線の周りに外積2回で作った座標系（線型代数 第5章）を使い、局所的な方向を$𝐝=x𝐭+y𝐛+z\widehat{𝐧}$とします。
23.9節のC++のパストレーサーも其の通りにしています。
シェーダーでは、23.7節と23.8節で使う2つの補助関数と一緒に書くと次の様になります：

static const float PI = 3.14159265;

// Cosine-weighted direction in the local frame (z = normal); pdf = cos(theta) / pi.
float3 sampleCosine(float2 u) {
  float r = sqrt(u.x);
  float phi = 2.0 * PI * u.y;
  return float3(r * cos(phi), r * sin(phi), sqrt(max(0.0, 1.0 - u.x)));
}

// Solid-angle pdf of a point picked uniformly on a light of area A (section 23.7).
float lightPdf(float dist, float cosLight, float area) {
  return dist * dist / (area * cosLight);
}

// Power heuristic (beta = 2): MIS weight of technique a against technique b (section 23.8).
float powerHeuristic(float pdfA, float pdfB) {
  float a = pdfA * pdfA;
  float b = pdfB * pdfB;
  return (a + b) > 0.0 ? a / (a + b) : 0.0;
}

const float PI = 3.14159265;

// Cosine-weighted direction in the local frame (z = normal); pdf = cos(theta) / pi.
vec3 sampleCosine(vec2 u) {
  float r = sqrt(u.x);
  float phi = 2.0 * PI * u.y;
  return vec3(r * cos(phi), r * sin(phi), sqrt(max(0.0, 1.0 - u.x)));
}

// Solid-angle pdf of a point picked uniformly on a light of area A (section 23.7).
float lightPdf(float dist, float cosLight, float area) {
  return dist * dist / (area * cosLight);
}

// Power heuristic (beta = 2): MIS weight of technique a against technique b (section 23.8).
float powerHeuristic(float pdfA, float pdfB) {
  float a = pdfA * pdfA;
  float b = pdfB * pdfB;
  return (a + b) > 0.0 ? a / (a + b) : 0.0;
}

static const float PI = 3.14159265;

// Cosine-weighted direction in the local frame (z = normal); pdf = cos(theta) / pi.
float3 sampleCosine(float2 u) {
  float r = sqrt(u.x);
  float phi = 2.0 * PI * u.y;
  return float3(r * cos(phi), r * sin(phi), sqrt(max(0.0, 1.0 - u.x)));
}

// Solid-angle pdf of a point picked uniformly on a light of area A (section 23.7).
float lightPdf(float dist, float cosLight, float area) {
  return dist * dist / (area * cosLight);
}

// Power heuristic (beta = 2): MIS weight of technique a against technique b (section 23.8).
float powerHeuristic(float pdfA, float pdfB) {
  float a = pdfA * pdfA;
  float b = pdfB * pdfB;
  return (a + b) > 0.0 ? a / (a + b) : 0.0;
}

23.5 重点的サンプリング

23.3節の分散の式をもう一度見て下さい：${\sigma }^2=\int f^2/p\mathrm{d}x-I^2$。
もし（正の$f$について）$p(x)=f(x)/I$と選べれば、全ての項が丁度$f/p=I$になり、分散は0です。
其れは出来ません。
探している答え$I$が必要だからです。
しかし方向は教えてくれます：pdfが被積分関数に似ている程、ノイズは少ない。
其の様に$p$を選ぶのが重点的サンプリングで、第12章は其の極端な例、$\cos \theta$のコサインサンプリングを見せました。

計算例：再び正弦

${\int }_0^{\pi }\sin x\mathrm{d}x$に対する3つの$p$の選び方：

三角形の分散${\int }_0^{\pi }{\sin }^{2}x/p(x)\mathrm{d}x-4=4.0670-4$（数値計算）は14分の1です：同じノイズに必要なサンプルは14分の1です。
同じ4つの乱数を三角形の逆CDFに通すと：

33%大き過ぎる代わりに9%小さいだけで、4つの項は全て2に近い値です。
理想のpdf $\frac{1}{2}\sin x$は逆関数法で$x=\arccos (1-2u)$とサンプリング出来、どのサンプルも丁度2を返します。

危険：小さ過ぎるpdf

重点的サンプリングは、事態を遥かに悪くする事もあります。
V字型のpdf $p(x)=4|x-\pi /2|/{\pi }^2$、つまり三角形を逆さにした物を取ります：正弦の頂上で丁度0です。
殆ど全ての所で正なので推定量は不偏の儘ですが、$\pi /2$の近くで項$\sin x/p(x)$が爆発し、$\int {\sin }^{2}x/p(x)\mathrm{d}x$は発散します：分散は無限大です。
${10}^7$個のサンプルの平均は其れでも2.000になりますが、標本分散（あたしの実行では32）はいつまでも落ち着かず、画像では稀に極端に明るいピクセルとして現れます。
規則：$f$が小さい所で$p$が大き過ぎるのは構わない（サンプルの無駄になるだけ）が、$f$が大きい所で$p$が小さ過ぎては絶対に駄目。

レンダラーでの重点的サンプリング

レンダリング方程式の被積分関数は積$f_r\cdot L_i\cdot \cos \theta$です。
入射光$L_i$は未知（再帰で計算する物）なので、分かっている因子で重点的サンプリングをします：

• BSDFサンプリング：出来るだけ$p\propto f_r\cos \theta$にする。ランバート面ならコサインサンプリングで、$f_r\cos \theta /p=\frac{\rho /\pi \cdot \cos \theta }{\cos \theta /\pi }=\rho$。
• 光源サンプリング：$p\propto$ ランプから直接届く光（23.7節）。
• 環境マップのサンプリング：HDRの空（第12章12.13節）の明るさの2次元の表を逆関数法で使う。行のCDF、次に選んだ行の中のピクセルのCDFで、太陽に殆どのサンプルが行く様にする。

第12章の一様な空からの放射照度なら、半球上の一様サンプリングの項は$2\pi L\cos \theta$で$\sigma =2\pi L/\sqrt{12}=1.814L$、コサインサンプリングでは毎回$\pi L$です。

GGXローブのサンプリング

光沢のあるマイクロファセットBRDF（第22章）は鏡面方向の周りに集中していて、コサインサンプリングは其処に殆どサンプルを送りません。
古典的な方法（Walter、Marschner、Li、Torrance、2007年）は、ハーフベクトル$𝐡$（マイクロファセットの法線）をGGX分布$D$からサンプリングし、視線方向を其の周りで反射させます。
$\alpha$をGGXの幅のパラメータとすると（多くのエンジンは$\alpha ={\text{roughness}}^2$とします）：

此のpdfは半球上で積分すると1になり（$D$は其の様に正規化されています）、CDFの逆関数から

次に${\omega }_i=2({\omega }_o\cdot 𝐡)𝐡-{\omega }_o$です。
$𝐡$の周りの反射は角度を2倍にするので、${\omega }_i$についてのpdfにはヤコビアンが付きます：

計算例。
$\alpha =0.3$、法線の真上から見て（${\omega }_o=𝐧$）、$u_1=0.5$：

（此処では${\omega }_o\cdot 𝐡=\cos {\theta }_h$なので$p({\omega }_i)=D/4$です。）
反射方向は法線から$2{\theta }_h=33.40\mathrm{°}$の方向に出て行きます。
ハーフベクトルの半分は法線から16.70°以内、90%は41.99°以内（$u_1=0.9$）に入ります：ローブはサンプリングされ、空っぽの半球はされません。
今日の改良版は、視点から実際に見える法線だけをサンプリングします（VNDFサンプリング、Eric Heitz、2018年）。
特に斜めから見る時に効きます。

見た目の違い

第12章の箱を、同じサンプル数（1ピクセルあたり16）で、拡散面での跳ね返りに3つのサンプリング方法を使ってパストレーシングした物です。
RMS誤差は1ピクセルあたり2048サンプルの参照画像に対して測りました：

コサインサンプリングは誤差を18%減らします：被積分関数の因子の1つに合っているからです。
大きく跳ぶのはランプも直接サンプリングした時（23.7節）で、最初の2枚のノイズの殆どは、小さなランプを運良く見つけた経路から来ているのです。

23.6 層化サンプリングと準モンテカルロ法

重点的サンプリングは、サンプルが平均してどこに行くかを変えます。
もう1つのレバーは、どれだけ均等に散らばるかです。
独立な乱数のサンプルは固まったり隙間を空けたりします。
23.3節の4つの内3つが山の頂上付近に集まった様にです。

層化サンプリング

層化サンプリングは、領域を$N$個の等しい部分（層）に分け、其々にランダムなサンプルを1つずつ置きます。
1次元では、サンプル$k$（0から数える）は$(k+u_k)/N$です。
同じ4つの数では：

$F_4=(1.4700+2.7410+2.9210+0.4671)/4=1.900$で、5%小さいだけです。
山の4分の1毎に丁度1つのサンプルが入るので、固まりようがありません。

推定値は不偏の儘で、分散は独立なサンプルより大きくなる事がありません。
滑らかな関数なら、誤差は1次元で$N^{-3/2}$、$d$次元で$N^{-1/2-1/d}$で減ります。
第2章（2.7節）の規則的な$4\times 4$のスーパーサンプリングの格子は、各サンプルをセルの中でランダムに動かすと層化サンプリングになります。
此れをジッターサンプリングと呼びます。

測ったRMS誤差（各点4000回）は予測通りです：1024サンプルで一様は0.030（$0.9669/32$）、三角形は0.0081（$0.2588/32$）、層化サンプリングは0.00006。
重点的サンプリングは線を下にずらし、層化は傾きを急にします。
両者は組み合わせられます：$u$の値を層化してから、逆CDFに通せば良いのです。

問題は又しても次元です。
$d$次元の完全なジッター格子には$k^d$個のサンプルが要ります：17次元で各軸に2層なら1ピクセルあたり131,072サンプルです。
其の為レンダラーは次元の小さな組（ピクセル、レンズ、跳ね返り毎）を別々に層化し、組同士の組み合わせを入れ替えます。
ラテン超方格やマルチジッターサンプリングの考え方です。

準モンテカルロ法

準モンテカルロ法（QMC）はランダムさを捨て、点が幾つでも空間を均等に埋める様に設計された決定的な低食い違い量列を使います。
点集合の食い違い量（ディスクレパンシー）は、其の最悪の偏りです：原点を角とする全ての箱について、中に入る点の割合と箱の体積の差の最大値です。
コクスマ・フラウカの不等式は誤差を（$f$の変動）×（食い違い量）で抑え、良い点列は約$(\log N{)}^d/N$、殆ど$1/N$に達します。

部品は基数$b$のラディカルインバース（radical inverse）です：番号$i$を$b$進法で書き、小数点を挟んで桁を鏡写しにします。
2進法では$1=1_2\to {0.1}_2=0.5$、$2={10}_2\to {0.01}_2=0.25$、$3={11}_2\to {0.11}_2=0.75$、$6={110}_2\to {0.011}_2=0.375$。
2進法の列$0.5,0.25,0.75,0.125,0.625,0.375,0.875,\dots$がファン・デル・コルプット列です：新しい点は、残っている一番大きな隙間の1つに入ります。

ハルトン列（John Halton、1960年）は次元毎に違う素数を基数にします：点$i$は（$i$の基数2のラディカルインバース、基数3の其れ、基数5の其れ、…）です。
計算例、2次元の点6：2進では上の通り0.375、$6={20}_3\to {0.02}_3=2/9=0.2222$なので、点は$(0.375,0.2222)$です。
3進の最初の8つは$1/3,2/3,1/9,4/9,7/9,2/9,5/9,8/9$です。

ソボル列（Ilya Sobol、1967年）は全ての次元で2進法を使い、次元毎に違う「方向数」をXORで組み合わせます。
ハルトン列（大きな基数では規則的な縞が出ます）より遥かに多くの次元で質を保ち、ビット演算で速いので、本格的なレンダラーでは普通こちらが選ばれます。

決定的な点列では誤差の見積もりが出来ません。
乱択準モンテカルロ法は、全ての点に同じランダムなずれを足して端で折り返したり（クランリー・パターソン回転）、桁をランダムに入れ替えたり（オーウェンのスクランブル）します。
すると各点は一様になって推定量は再び不偏になり、点は均等な広がりを保ちます。
滑らかな2次元の被積分関数$f(x,y)=e^{-(x+y)}$で、2000回の乱択について測ったRMS誤差は：

1024サンプルでは、ずらしたハルトン列は乱数より13倍正確です。
此処でジッター格子が勝つのは、2次元で$N$が平方数だからに過ぎません。
ハルトン列はどんな$N$でも使え、1サンプルずつ追加出来、多くの次元に広げられます。
1ピクセルあたり1、2サンプルしか取らないリアルタイムのレンダラーでは、残る誤差がどう見えるかの方が大事です：ブルーノイズのパターンは誤差を細かいざらつきに変え、目にもデノイザーにも塊より扱い易くします。

23.7 光源のサンプリング

コサインサンプリングでも何故あれ程ノイズが出たのでしょうか？
ランプが小さく、コサインでサンプリングした方向が滅多に当たらないからです。

ランダムな反射は何回ランプを見つけるか

0.6 m × 0.6 mの四角いランプ（$A=0.36$ m²、第12章の箱のランプ）、放射輝度$L_e=14$の中心の、2 m真下の床の点を考えます。
数値積分すると、ランプからの放射照度は$E=1.2234$ W/m²です。
コサインでサンプリングした方向がランプに当たる確率は

ランプを見るのは跳ね返りの2.8%だけで、推定値は0か大きな値かのどちらかになります。
最悪の種類のノイズです。
点光源（立体角0）なら、確率は丁度0です。

光源上の点をサンプリングする

光源サンプリングは此れを逆にします：ランプ上の点$𝐲$を面積について一様に選び（$p_A=1/A$）、其処へシャドウレイを送ります。
此れが次イベント推定（NEE）で、第12章のパストレーサーも既にしていました。
方向とランプの面積は、次の図の幾何で結ばれます：

第12章の立体角と同じく、パッチは遠い程、又傾いている程小さく見えます。
従って面積についてのpdf $1/A$は、$𝐱$での立体角については

で、1サンプルの直接光の推定量は

$V$はシャドウレイが遮られなければ1、遮られれば0です。

計算例。
ランプの下の床の点で、床は$\rho =0.75$の白なので$f_r=0.75/\pi$です。
ランプの中心のサンプル：$r=2$、どちらのコサインも1で、放射照度の部分は$14\times 1\times 1\times 0.36/4=1.26$。
床の点から見て$(0.3,2,0.3)$の角のサンプル：$r^2=0.09+4+0.09=4.18$、$\cos \theta =\cos {\theta }_l=2/\sqrt{4.18}=0.9782$で、放射照度の部分は$14\times 0.9569\times 0.36/4.18=1.1538$。
どのサンプルも此の2つの値の間に入り、真の値1.2234に近くなります。
コサインサンプリングの「0か、36倍大き過ぎるか」とは大違いです。
中心のサンプルからの反射放射輝度は$\frac{0.75}{\pi }\times 1.26=0.3008$です（真の値は$\frac{0.75}{\pi }\times 1.2234=0.2921$）。
立体角で言えば、中心方向の光源のpdfは1ステラジアンあたり$p_{\omega }=4/0.36=11.11$、コサインサンプリングは其処で$1/\pi =0.318$です：光のある所に35倍多くのサンプルが行きます。

沢山の光源

ランプが何千個もあるなら、確率$P_j$で1つをランダムに選び、其の寄与を$P_j$で割ります。
1段上の$f/p$の規則です。
電力に比例させると、100 W、50 W、10 Wのランプは$P=0.625,0.3125,0.0625$で、弱いランプが選ばれた時は其の寄与が16倍されます。
本格的なレンダラーは明るくて近い光源を優先する木構造（「ライトBVH」）から光源を選び、ReSTIR（第24章）は良い光源サンプルをピクセル間やフレーム間で使い回し、何百万もの光源をリアルタイムで扱います。

23.8 多重重点的サンプリング

光源サンプリングは、小さなランプの下の艶消しの床には最適です。
では床を殆ど鏡にし、ランプを巨大な窓にしてみましょう。
光源のサンプルは窓全体に散らばり、殆どは鋭いBRDFのローブの外（値は殆ど0）、稀に1つが中（巨大な値）に入ります。
一方BSDFサンプリングは全てのレイをローブの中に送り、其処では大きな窓に簡単に当たります。
どちらの手法が勝つかは材質と光源次第で、本物のシーンには全ての組み合わせが同時にあります。

2つの手法を組み合わせる

多重重点的サンプリング（MIS）は1995年にEric VeachとLeonidas Guibasが発表した物で（"Optimally Combining Sampling Techniques for Monte Carlo Rendering"）、両方の手法を使い、各サンプルを其の点で其の手法がどれだけ得意かで重み付けします。
手法$i$からpdf $p_i$で$n_i$個のサンプル$X_{i,j}$を取ると：

$f(x)\ne 0$となる全ての$x$で重みの和が1、${\sum }_i{w}_i(x)=1$なら（且つ$p_i=0$の所では$w_i=0$なら）、此れは不偏です。
証明は1行です：期待値は${\sum }_i\int w_{i}f\mathrm{d}x=\int ({\sum }_i{w}_i)f\mathrm{d}x=\int f\mathrm{d}x$。
問題は重みの選び方です。
VeachとGuibasは2つのヒューリスティックを提案しました：

普通は$\beta =2$です。
其の点でpdfが高い手法（そういう点をよく見つける手法）は大きな重みを貰い、運だけで見つかったサンプル（pdfが低く、$f/p$が巨大）は小さな重みを貰います。
其れこそが、放っておけばファイアフライになるサンプルです。
バランスヒューリスティックでは各項が$f(x)/{\sum }_k{n}_k{p}_k(x)$になり、全てのサンプルが1つの混合pdfから来た様に見えます。
Veachは、其の分散が他のどんな重み付けよりも高々$(1/{min}_i{n}_i-1/{\sum }_i{n}_i)I^2$しか大きくならない事を証明しました。
パワーヒューリスティックは重みを鋭くし、1つの手法が明らかに正しい時に効きます。
多くのレンダラーが使うのは此方です。

計算例：重み

23.7節の艶消しの床で、ランプの中心への方向、各手法から1サンプルずつ：$p_l=11.11$、$p_b=0.3183$。

此処では明らかに光源サンプリングが正しい手法で、MISは殆ど全ての重みを其方に渡します。

光沢のある面で、同じ様なランプについて、或る方向で$p_l=2$、$p_b=8$の場合：

どちらの場合も重みの組は足して1になり、其れが推定量を不偏に保ちます。

計算例：分散

1次元の「フラットランド」版で、効果を測る事が出来ます。
20°の周りの幅$s$の釣鐘型の光沢ローブ（BSDFサンプリングが正確に従う）を、角度の或る区間を覆うランプ（光源サンプリングが一様にサンプリングする）が照らし、被積分関数はランプの中でローブ×コサインです。
全部で2サンプル（1つの手法で2つ、又はMISで1 + 1）の時の、正確な相対標準偏差（標準偏差 / 真の値）：

単独の手法は、一方の場合では非常に良く、もう一方では酷い物です。
MISは決して酷くなりません：正しい手法（2サンプルの代わりに1サンプル）には負けますが、150〜180%という大惨事を無くし、どちらの場合もパワーヒューリスティックがバランスヒューリスティックに勝ちます。
レンダラーは各ピクセルがどちらの場合なのかを知らないので、此の頑健さこそが必要なのです。

中間の場合（0°の周りの幅12°のローブ、−5°から35°のランプ）の2つのpdfと重みです：

ローブが高い所ではBSDFサンプリングが重みを取り、ローブの縁に向かって光源サンプリングが引き継ぎます。
重みが交わるのは9°、丁度$p_b=p_l$となる所です。
MISには1サンプル版もあります：確率$c_i$で手法$i$を選び、1サンプル取って$w_{i}f/(c_i{p}_i)$を使います。
材質の層（拡散のベースかクリアコートか）を選ぶ時に便利です。

23.9 パストレーシングの詳細

パストレーサーが何を推定しているのかを考え、23.10節の手法を組み立てる為に、Veachの博士論文（"Robust Monte Carlo Methods for Light Transport Simulation"、スタンフォード大学、1997年）はレンダリング方程式を経路についての普通の積分に書き直しました。

経路積分

長さ$k$の経路は点の並び$\overline{x}=({𝐱}_0,{𝐱}_1,\dots ,{𝐱}_k)$です：${𝐱}_0$はカメラ上、${𝐱}_k$は光源上、他は面の上にあります。
ピクセル$j$の値は

全ての経路の長さについての、全ての頂点の位置についての積分です。
測定寄与関数$f_j$は、放射される光、跳ね返る毎のBSDF、カメラの感度$W_e$、そして隣り合う頂点の間の幾何項を掛け合わせた物です：

$G$は23.7節の$\cos {\theta }_l/r^2$の面と面の版です：2つの点がお互いをどれだけ「見ている」か。
2 m離れ、コサインが0.8と0.6で遮る物の無い2点では$G=0.8\times 0.6/4=0.12$（1 m²あたり）です。
（第12章のラジオシティ法の形態係数は、此の$G$を$\pi$で割って2つのパッチで平均した物です。）
ヘックバートの表記にも正確な意味が付きます（LDSEは頂点が光源、拡散面、鏡面、目の経路の集合）。
そして全てのレンダリングアルゴリズムは、経路$\overline{x}$をサンプリングし、$f_j(\overline{x})$を其れを作った確率密度で割る方法の違いなのです。

スループット

パストレーサーは各頂点で点ではなく方向を選びます。
立体角のpdf $p_{\omega }$の方向は、次の頂点に面積のpdf $p_{\omega }\cos \theta \prime /r^2$を与え（又しても23.7節）、其れが$G$の同じ因子と打ち消し合って、スループットの更新式が残ります：

$\beta$はカメラで1から始まり、経路が放射される光$L_e$に届く度に、ピクセルは$\beta L_e$を受け取ります。
コサインサンプリングするランバート面では因子はアルベド$\rho$なので（23.5節）、$k$回の拡散反射の後は$\beta ={\rho }_1{\rho }_2\cdots {\rho }_k$です。

ロシアンルーレットと其の代償

経路を決まった最大の深さで止めるとバイアスが出ます：もっと跳ね返る必要のある光が失われるのです。
ロシアンルーレット（1990年にJames ArvoとDavid Kirkがグラフィックスに持ち込みました）は確率$q$で続けて$\beta$を$q$で割り、其れ以外では止めます。
残りの寄与$C$は、確率$q$で$C/q$に、其れ以外では0になります：

不偏ですが、只ではありません。
増える分散は

第12章の例（$C=0.5$、$q=0.7$）では$0.25\times 0.4286=0.107$で、其の代わりに仕事の30%を省けます。
大事なのは効率$1/(\text{分散}\times \text{時間})$です：残りの寄与が小さい経路を止める時、ロシアンルーレットは得になります。
其の為、普通は何回か跳ね返った後で$q=min(1,max({\beta }_R,{\beta }_G,{\beta }_B))$とします：暗い経路は間引かれ、生き残った経路はスループットが1近く迄持ち上げられます。
一定の$q$なら、経路は平均$1/(1-q)$回跳ね返ります（$q=0.8$なら5回）。

小さなパストレーサー

答えの分かっているシーンでの、C++のパストレーサーの中心のループです：カメラは閉じた球の中にあり、球の内壁はアルベド$\rho$のランバート面で、放射輝度$L_e$で光っています。
第12章の等比級数（12.7節）から、正確な答えは$L_e(1+\rho +{\rho }^2+\dots )=L_e/(1-\rho )$、$\rho =0.8$なら5です。
ループは球との交差を求め、放射にスループットを掛けて足し、コサイン重み付きの方向をサンプリングし、スループットに$\rho$を掛け、ロシアンルーレットをします。
線型代数シリーズのlinalg.hhのVec3、dot、cross、normalizeを使います：

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <random>
#include "linalg.hh"

// A camera at the center of a closed sphere (radius 1). The inside wall is
// Lambertian with albedo rho and emits radiance Le. Exact answer: Le / (1 - rho).
const float PI = 3.14159265f;
std::mt19937 gen(7);
float rnd() { return (gen() >> 8) * (1.0f / 16777216.0f); }

// Cosine-weighted direction around the unit normal n (pdf = cos / pi).
Vec3 sampleCosine(Vec3 n) {
  float u1 = rnd(), u2 = rnd();
  float r = std::sqrt(u1), phi = 2.0f * PI * u2;
  Vec3 a = std::fabs(n.x) > 0.9f ? Vec3{0, 1, 0} : Vec3{1, 0, 0};
  Vec3 t = normalize(cross(a, n));
  Vec3 b = cross(n, t);
  return r * std::cos(phi) * t + r * std::sin(phi) * b + std::sqrt(1.0f - u1) * n;
}

// The core loop of a path tracer. maxDepth < 0 means: no limit, Russian roulette instead.
float radiance(Vec3 o, Vec3 d, float rho, float Le, int maxDepth, int *bounces) {
  float L = 0.0f, beta = 1.0f;  // beta = path throughput
  for (int depth = 0;; depth++) {
    float b = dot(o, d), c = dot(o, o) - 1.0f;  // hit the sphere from inside
    Vec3 p = o + (-b + std::sqrt(b * b - c)) * d;
    Vec3 n = -1.0f * p;                         // inward unit normal
    L += beta * Le;                             // emission at the hit
    *bounces = depth + 1;
    if (maxDepth >= 0 && depth >= maxDepth) break;
    d = sampleCosine(n);
    beta *= rho;                                // f cos / p = (rho/pi) cos / (cos/pi)
    if (maxDepth < 0) {                         // Russian roulette
      float q = std::min(beta, 1.0f);
      if (rnd() >= q) break;
      beta /= q;
    }
    o = p + 1e-4f * n;
  }
  return L;
}

void run(const char *name, float rho, int maxDepth) {
  const int paths = 1000000;
  double sum = 0, sq = 0, steps = 0;
  for (int i = 0; i < paths; i++) {
    int k = 0;
    double v = radiance({0, 0, 0}, {0, 0, -1}, rho, 1.0f, maxDepth, &k);
    sum += v;
    sq += v * v;
    steps += k;
  }
  double mean = sum / paths, sd = std::sqrt(std::max(0.0, sq / paths - mean * mean));
  printf("%-22s mean = %.4f  sd per path = %.3f  hits per path = %.3f\n", name, mean, sd, steps / paths);
}

int main() {
  float rho = 0.8f;
  printf("exact Le / (1 - rho)   = %.4f\n", 1.0f / (1.0f - rho));
  run("max depth 4", rho, 4);
  run("max depth 16", rho, 16);
  run("Russian roulette", rho, -1);
  return 0;
}

exact Le / (1 - rho)   = 5.0000
max depth 4            mean = 3.3616  sd per path = 0.000  hits per path = 5.000
max depth 16           mean = 4.8874  sd per path = 0.000  hits per path = 17.000
Russian roulette       mean = 5.0006  sd per path = 4.478  hits per path = 5.001

3行でバイアスと分散が分かります。
5回当たった所で止めると、此の対称なシーンではノイズは0ですが、33%暗過ぎます（$1+0.8+\dots +{0.8}^4=3.3616$）。
17回でも2%暗過ぎます（$(1-{0.8}^{17})/0.2=4.8874$）。
ロシアンルーレットは平均で正しく（5.0006、100万本の経路の標準誤差は$4.478/1000=0.0045$）、1本あたり5回、$1/(1-0.8)$回しか当たりません。
しかし1本1本の経路はノイズを持ちます：其の値は当たった回数で、標準偏差$\sqrt{0.8}/0.2=4.472$の幾何分布の確率変数です。
此れが第12章（12.13節）のホワイトファーネステストで、真面目なレンダラーは皆此れを走らせます。

鏡とガラス：デルタ分布

完全な鏡のBRDFはデルタ分布です：1つの方向以外では0で、合計は$F$（フレネル反射率）。
パストレーサーは鏡面方向を取り、$f_s\cos \theta /p$を其の極限$F$に置き換えます（第12章のレンダラーの様に確率$F$と$1-F$で反射か屈折かを選ぶ時は1に）。
其処では光源サンプリングは役に立たない（ランプ上のランダムな点が丁度鏡面方向にある確率は0）ので、鏡面の頂点ではNEEを飛ばし、鏡面反射の直後にランプに当たった時は、其の放射をMISの重み1で数えます。
此の章のレンダラーでは、フラグ1つです。
此れで第12章のガラス球の下の頑固な斑点の理由が分かります：目から辿るコースティクスの経路LSDEは、床で跳ね返ってガラスを通り、小さなランプに偶然当たらなければならず、点光源なら其れは絶対に起きないのです。
LS+DS+Eの経路（鏡に映ったコースティクス、水面越しに見るプールの底）は、拡散面の頂点で繋ぐ全ての手法にとって最も難しい物です。

パストレーサーの中のMIS

23.7節と23.8節を合わせると、鏡面でない全ての頂点は2つの事をします：

at a diffuse or glossy vertex x (throughput beta, outgoing direction wo):
    light sample:  pick y on a light, pdf pL (solid angle)
                   if visible: add beta * f(wo, wi) * cos * Le / pL * w(pL, pB(wi))
    BSDF sample:   pick wi with pdf pB, beta *= f(wo, wi) * cos / pB
                   trace the ray; if it hits a light:
                       add beta * Le * w(pB, pL(hit))
                   otherwise continue the path from the hit

w(a, b)は23.4節のシェーダーのパワーヒューリスティックです。
此れが23.5節の図の「コサイン＋光源（MIS）」のレンダラーです。

23.10 双方向の手法、メトロポリス、フォトンマッピング

MIS付きのパストレーシングは映画制作の標準ですが、小さな光源からのコースティクス、隙間から入る光、ガラス越しに見る光は難しい儘です。
主な代替手法の考え方を見ます。

双方向パストレーシング

双方向パストレーシング（BDPT。Eric LafortuneとYves Willemsが1993年に、独立にVeachとGuibasが1994年に発表）は、カメラから1本、光源から1本の部分経路を辿り、両者の頂点をシャドウレイで繋ぎます：

同じ完全な経路は、幾つもの方法で作れます：光源側から$s$個、視点側から$t$個の頂点を取ると、$n$頂点の経路は$n+1$通り（$s+t=n$で$s=0,1,\dots ,n$）に作れます。
図の5頂点（${𝐱}_0,{𝐱}_1,{𝐱}_2,{𝐲}_1,{𝐲}_0$）なら6つの手法です：$s=0$は純粋なパストレーシング（視点側の経路が自分でランプに当たる）、$s=1$は次イベント推定付きのパストレーシング、$s=2$は図の接続、と続き、$s=5$は光源側の経路が自分でカメラに当たる物です（ピンホールカメラでは不可能なので、pdfは0）。
全ての$n+1$個のpdfに対するMISが其れらを組み合わせます。
MISは正に此の為に考えられたのです。
光源側の頂点を直接カメラに繋ぐ$t=1$の手法がライトトレーシングで、直接見えるコースティクス（LS+DE）を簡単に描けます。
代償は複雑さと、1サンプルあたりの仕事の多さです。

メトロポリス光輸送

BDPTにとっても稀な経路があります。
半開きのドアから入る光の様な物です。
メトロポリス光輸送（MLT。VeachとGuibas、1997年）はメトロポリス・ヘイスティングス法（1953年、1970年）を使います：光を運ぶ経路が一度見つかったら、其の近くを探索するのです。
経路のマルコフ連鎖を作ります：今の経路$\overline{x}$から変異$\overline{y}$（頂点を動かす、部分経路を取り替える…）を提案し、次の確率で受け入れます。

$f$は経路の明るさ、$T$は一方から他方を提案する確率です。
対称な変異なら$T$は打ち消し合います：明るさ2の経路から0.5の経路へは確率$0.5/2=25\%$で、3の経路へは必ず移ります。
連鎖は経路を其の明るさに比例して訪れ、画像全体の明るさは別の普通の推定で決めます。
連続するサンプルは相関するので、ノイズは染みの様になり、収束にむらが出ます。
又、連鎖の出発点が最初の内は結果を偏らせます。
主標本空間MLT（Csaba Kelemenら、2002年）は、普通のパストレーサーを動かす乱数の方を変異させる事で、此の手法を大幅に簡単にしました。

プログレッシブフォトンマッピング

第12章のフォトンマッピングは、フォトンの密度$\sum {\mathrm{\Phi }}_p/(\pi r^2)$から放射輝度を推定します。
其れは半径$r$の範囲でぼけるのでバイアスがあり、フォトン数が無限大に、同時に$r$が0に向かう時にだけ一致性を持ちます。
1つのフォトンマップを記憶しておく方式では、其れは出来ません。
プログレッシブフォトンマッピング（Toshiya Hachisuka、Shinji Ogaki、Henrik Wann Jensen、2008年）は、フォトンを記憶せずに何回でもフォトンのパスを走らせます。
目からのレイの各交点は半径$R_i$と数$N_i$を持ちます。
或るパスで半径の中に新しいフォトンが$M_i$個見つかったら、其の割合$\alpha$だけを残し、密度が変わらない様に円盤を縮めます：

論文は$\alpha =0.7$を使っています。
計算例：$N_i=100$、$M_i=50$なら$N_{i+1}=135$、$R_{i+1}=R_i\sqrt{135/150}=0.9487R_i$です。
半径はパス毎に縮み、バイアスは0に向かい、パストレーシングを打ち負かすSDSの経路も含めて画像は収束します。
確率的プログレッシブフォトンマッピング（HachisukaとJensen、2009年）は其れを光沢反射や被写界深度に広げ、2012年には2つのグループがフォトンの密度推定とBDPTをMISの下で組み合わせました：vertex connection and merging（VCM。Iliyan Georgievら）とunified path sampling（UPS。Hachisukaら）です。
ボリュームでは媒質の中の散乱の頂点が加わり、23.4節の自由行程のサンプリングを使います（第28章）。

23.11 バイアス・分散・一致性

此処迄使ってきた3つの言葉を、はっきりさせておきます。

• バイアス（偏り）：$E[F_N]-I$、平均しても残る誤差。バイアスが0の推定量が不偏。
• 分散：$E[F_N]$の周りのランダムな散らばりで、ノイズとして見える。
• 一致性：仕事（サンプル、フォトン、パス）を無限に増やすと誤差が0に向かう事。分散が有限の不偏推定量には一致性があり、バイアスのある推定量も、バイアスが仕事と共に縮めば一致性を持てる。

実際に目に見えるのは全体の誤差で、平均二乗誤差（MSE）で測ります：

計算例。
同じレンダリング時間で、レンダラーAは不偏で分散0.04、レンダラーBはバイアス0.05で分散0.01とします。
A：MSE = 0.04。
B：MSE = $0.01+{0.05}^2=0.0125$で、約3倍良い。
16倍長くレンダリングすると、AのMSEは$0.04/16=0.0025$に下がりますが、Bはバイアスが残る（Bに一致性が無ければ）ので$0.01/16+0.0025=0.0031$で、其処からはAが勝ちます。
バイアスのある手法（フォトンマッピング、クランプ、デノイズ）は早く綺麗に見え、不偏な手法は極限で勝ち、誤差をサンプルから見積もれるので信頼し易いのです。

ファイアフライとクランプ

ファイアフライとは、周りよりずっと明るい1つのピクセルです：pdfが小さく$f$が大きいサンプル1つで、典型的には目から見つけたコースティクスです。
サンプルの値が確率99.9%で0.5、0.1%で500のピクセルを考えます：

真の明るさの半分は、稀な経路から来ています。
1サンプルは$E[X^2]=0.999\times 0.25+0.001\times 250000=250.25$、分散$250.25-{0.9995}^2=249.25$、標準偏差15.8で、1ピクセルあたり64サンプルでも$15.8/8=1.97$、ピクセルの値の2倍です。
見える物：此の様なピクセルの$1-{0.999}^{64}=6.2\%$は稀なサンプルを捕まえて約$(63\times 0.5+500)/64=8.3$を示し、白い点になります。
残りの93.8%は0.5を示し、本来の半分の明るさです。

クランプは全てのサンプルを最大値、例えば10に制限します。
ファイアフライは消えますが、平均はサンプル数に関係なく0.9995ではなく$0.4995+0.001\times 10=0.5095$になります。
本格的なレンダラーも其れでもクランプします（多くは間接光だけ。下のデノイズの図も其うしています）。
コースティクスの明るさが欠ける方が、斑点より見る人は気にならないからです。
本当の解決は、其の経路をよく見つける手法です（23.10節）。

23.12 デノイズと適応的サンプリング

映画のフレームは今でも1ピクセルあたり数百から数千サンプルを必要とし、ゲームのフレームに使えるのは1サンプル程度です。
デノイズはノイズのある画像をフィルタで処理して此の差を埋めます。
エッジやテクスチャをぼかさない為に、レンダラーがただで持っている情報を使います。

Gバッファを使う

レンダラーは各ピクセルの最初の交点の法線、位置又は深度、アルベドを知っています。
ディファードレンダラーのGバッファです。
此れらは殆どノイズが無く、エッジがどこにあるかを正確に示します。
クロスバイラテラルフィルタは、第15章のバイラテラルフィルタのエッジの重みを別の画像から取る物で、法線、位置、色が似た近傍だけと平均します。
更に2つの工夫があります：

• アルベドの分離：フィルタの前に色をアルベドで割り、後で掛け戻す。フィルタは照明を均し、テクスチャはくっきりした儘になる。
• 光源を外す：直接見えるランプにはノイズが無く、周りに滲ませてはいけない。

エッジ回避à-trousフィルタ

サンプルが少ないと幅の広いフィルタが要りますが、65 × 65のクロスバイラテラルの窓は1ピクセルあたり4,225タップかかります。
エッジ回避à-trousウェーブレット変換（Holger Dammertz、Daniel Sewtz、Johannes Hanika、Hendrik Lensch、2010年）は、5 × 5のフィルタを、毎回タップの間隔を2倍に広げながら何回も掛けます（「à trous」はフランス語で「穴の開いた」）：

• 1次元の重みはB3スプラインの$\frac{1}{16},\frac{1}{4},\frac{3}{8},\frac{1}{4},\frac{1}{16}$で、5 × 5のカーネルは其の積、25タップ。
• パス$i$（0から）は$2^i$ピクセル離れたタップを使うので、5回目のパスは片側$2\times 16=32$ピクセル先に届き、65 × 65の窓になる。
• 各タップの重みには色、法線、位置の差からのエッジ停止関数が掛かり、フィルタはエッジを越えない。

5回のパスで1ピクセルあたり4,225タップではなく$5\times 25=125$タップです。
論文は、残るノイズが小さくなる分、パス毎に色の許容幅を半分にしています。

計算例：1つのタップの重み。
2回目のパスでカーネル座標(+1, 0)、つまり2ピクセル右の近傍のスプラインの重みは$\frac{1}{4}\times \frac{3}{8}=0.09375$です。
SVGFのエッジ停止関数（下）は、此れに法線の重み$(𝐧\cdot 𝐧\prime {)}^{128}$と輝度の重み$\exp (-|l-l\prime |/(4s))$を掛けます。
$s$はノイズの局所的な標準偏差です。
$𝐧\cdot 𝐧\prime =0.99$なら法線の重みは${0.99}^{128}=0.276$、$𝐧\cdot 𝐧\prime =0.9$の折れ目を越えると${0.9}^{128}=1.4\times {10}^{-6}$で、実質0です。
$s=0.05$で輝度の差が0.1なら$\exp (-0.1/0.2)=0.607$で、タップの重みは正規化前で$0.09375\times 0.276\times 0.607=0.0157$です。

SVGF：時間と分散を加える

SVGF（spatiotemporal variance-guided filtering。Christoph SchiedらNVIDIAとKITのグループ、2017年）は、此れを1ピクセルあたり1サンプルのリアルタイムのデノイザーにしました：

1. 時間方向の蓄積：前のフレームの結果をモーションベクトルで再投影して混ぜる。$c=(1-\alpha )c_{\text{history}}+\alpha c_{\text{new}}$、$\alpha =0.2$。此の移動平均はノイズの分散を$\alpha /(2-\alpha )=0.111$に減らし、9フレームの平均と同じ位。遮蔽が外れた所では履歴をリセットする。
2. 分散の推定：時間方向に蓄積したモーメントから、ピクセル毎の輝度の分散を求める（履歴が短い所では空間的に推定）。
3. à-trousフィルタ：法線、深度、分散で導かれる輝度の重みを使って5回。ノイズの多い所は強く均し、収束した所には殆ど触れない。分散も色と一緒にフィルタする。

アルベドを分離して、論文は当時のGPUで1920 × 1080に約10 msと報告しています。

1ピクセルあたり8サンプル（間接光は2でクランプ）の箱のシーンを、ガウスぼかしと、此の考え方のà-trousフィルタ（5回のパス。法線、位置、材質によるエッジ停止。空間的な分散の推定で導かれる輝度の重み。アルベドの分離。ランプは除外）で処理した物です。
RMS誤差は表示されるsRGBの値で測りました：

ガウスぼかしは入力より悪くなります：全ての所でノイズをぼけと交換し、ランプを滲ませるのです。
エッジを考えるフィルタは誤差を3分の1減らし、形をくっきり保ちます。
しかしガラス球の下のコースティクスは消えています。
クランプされ、平均されて無くなったのです：デノイズのバイアスが目に見える形です。

ニューラルデノイザー

2017年以降、最良のデノイザーは、ノイズのあるレンダリングと収束したレンダリングの組で学習し、Gバッファも入力に取るニューラルネットワークです。
例えばChakravarty Alla Chaitanyaらの再帰型オートエンコーダーや、Steve Bako、Thijs Vogelsらのカーネル予測ネットワーク（どちらも2017年）です。
執筆時点では、IntelのOpen Image DenoiseとNVIDIAのOptiXのデノイザーがレンダラーで広く使われています。
ゲームではNVIDIAのDLSS Ray Reconstruction（2023年）が手作業で調整したデノイザーを置き換え、デノイズとアップスケーリングを一緒に行います。
ニューラルレンダリングは第29章の話題です。

適応的サンプリング

照らされた壁は数サンプルで収束しますが、コースティクスの光沢のある映り込みには何千サンプルも要ります。
適応的サンプリングは各ピクセルの平均$\mu$と標本標準偏差$s$を見積もり、相対誤差$s/(\sqrt{N}\mu )$が目標$t$を下回る迄サンプリングを続けます：

計算例。
目標$t=2\%$。
$\mu =0.3$、$s=0.6$のノイズの多いピクセルには$N\ge (0.6/(0.02\times 0.3){)}^2={100}^2=10,000$サンプル、$s=0.1$の落ち着いたピクセルには$(0.1/0.006{)}^2=278$サンプルで足ります。
分散の推定自体にもノイズがあり（稀なコースティクスに未だ出会っていないピクセルは収束した様に見えるので、最小のサンプル数を決めておきます）、ピクセル自身のサンプルで止める事は僅かなバイアスを生むので、頑健なレンダラーは別のサンプルバッファから誤差を見積もります。
レンダラーは此れを「ノイズの閾値」として見せています。
リアルタイムのレイトレーシングは、更にサンプルをピクセル間やフレーム間で使い回します。
ReSTIRがしている事です（第24章）。

23.13 まとめ

• レンダリングの積分は次元が多い。格子による公式は其処で破綻するが、モンテカルロ法の誤差$\sigma /\sqrt{N}$は次元に依存しない。独立な分散は足し算で、定数倍は2乗で効くからである。
• 推定量$\frac{1}{N}\sum f(X_k)/p(X_k)$は、$f\ne 0$の所で$p>0$なら不偏。1サンプルあたりの分散は$\int f^2/p-I^2$。
• サンプリング：逆関数法$X=P^{-1}(U)$（指数分布の自由行程$-\ln (1-u)/\sigma$、コサイン半球$\sin \theta =\sqrt{u_1}$）、代わりの手段として棄却法。
• 重点的サンプリング：$f$に似た形の$p$は分散を下げる（正弦なら三角形で14分の1）。$f$が大きい所で$p$が小さ過ぎると無限大にもなる。レンダラーはBSDF（コサイン、$p({\omega }_i)=D\cos {\theta }_h/(4{\omega }_o\cdot 𝐡)$のGGXのハーフベクトル）、光源、環境マップをサンプリングする。
• 層化（ジッター）サンプリングと準モンテカルロ法（ハルトン列1960年、ソボル列1967年）はサンプルを均等に広げる。乱択QMCは不偏性を保つ。
• 光源サンプリング（NEE）：$p_{\omega }=r^2/(A\cos {\theta }_l)$。点光源に届く唯一の方法。光源が多ければ確率で1つを選び、其の確率で割る。
• MISは和が1になる重みで手法を組み合わせる：バランスヒューリスティック$p_i/\sum p_k$、2乗を使うパワーヒューリスティック。決して破綻せず、最良の手法に近い。
• 経路積分：ピクセル = 経路についての$L_e\cdot G\cdot f_s\cdots W_e$の積分。パストレーサーは跳ね返る毎にスループットに$f\cos /p$を掛ける。ロシアンルーレットは不偏だが分散$C^2(1/q-1)$を加え、固定の深さの上限にはバイアスがある。鏡とガラスはデルタBSDFで、其処では光源サンプリングをしない。
• BDPTは視点側と光源側の部分経路を繋ぐ（$n$頂点に$n+1$通りの手法、MISで組み合わせる）。MLTは明るい経路を変異させる。プログレッシブフォトンマッピングは半径を縮め、バイアスはあるが一致性がある。
• MSE = 分散 + バイアス²。クランプはバイアスと引き換えにファイアフライを消す。デノイザーはGバッファを使う（クロスバイラテラル、à-trous、SVGF、ニューラルネットワーク）。適応的サンプリングは見積もった誤差が大きい所にサンプルを使う。

23.14 用語

23.15 課題

1. ★ 一様なサンプル$u=0.2,0.5,0.9$で${\int }_0^1{x}^2\mathrm{d}x$を推定し、次に同じ$u$から逆関数法でpdf $p(x)=2x$に従うサンプルを作って推定して下さい。両方の推定量の1サンプルあたりの正確な分散も求めて下さい。（確認：正確な積分は$1/3$、2つ目の推定量の分散は$1/72$。）
2. ★ 基数2と3でハルトン列の点11を、又ファン・デル・コルプット列の点9から12を求めて下さい。2進の最初の12個の値を数直線に記し、幅1/8の8つの区間の内、8点の後と12点の後に其々幾つの区間に点があるか数えて下さい。
3. ★★ $\alpha =0.5$、$u_1=0.25$のGGXサンプリングで、${\theta }_h$、$D(𝐡)$、$p(𝐡)$、そして法線の真上から見る時の$p({\omega }_i)$と反射方向の角度を求めて下さい。$u_1=0.75$でも計算して比べて下さい。
4. ★★ 光源のサンプルが$p_l=6$で、同じ方向のBSDFのpdfが$p_b=2$です。BSDFのサンプルは$p_b=3$で、$p_l=0.5$の所で光源に当たります。全てのバランスとパワーの重みを求め、各点で和が1になる事を確かめて下さい。次に、スループット0.2の経路が$q=0.2$のロシアンルーレットに出会います。取り得るスループットを挙げ、期待値が変わらない事を示し、残りの寄与$C=0.2$について増える分散を求めて下さい。
5. ★★★ 点光源がガラス球を通して艶消しのテーブルを照らし、カメラは其のコースティクスを鏡の中にしか見ません。経路をヘックバートの表記で書き、NEE付きのパストレーシング、ライトトレーシング、BDPT、プログレッシブフォトンマッピングのどれが其れを描けるか、其の理由と共に説明して下さい。

選択問題

問1 ★ 次のうち、モンテカルロ推定量$F_N=\frac{1}{N}\sum f(X_k)/p(X_k)$について正しい物はどれか。

問2 ★★ 次のうち、多重重点的サンプリングについて正しい物はどれか。

問3 ★★ 次のうち、バイアスと一致性について正しい物はどれか。

23.16 実践課題

1つの点の直接光を、幾つかのサンプリング方法で手計算します。
方眼紙、定規、分度器、電卓を用意して下さい。
フラットランド、つまり2次元の断面で考えます：方向は法線から測った角度$\theta$（右が正）で、「放射照度」は$-90\mathrm{°}$から$90\mathrm{°}$迄の$E=\int L(\theta )\cos \theta \mathrm{d}\theta$です。
立体角の代わりに、ラジアンの$\mathrm{d}\theta$を使います。
点は床の原点（法線は上向き）、ランプは高さ$h=2$で$x=0$から$x=2$迄（長さ$\mathrm{\ell }=2$）の水平な線分で、下を向き、放射輝度は$L=1$です。
どの手法にも乱数$u=0.10,0.35,0.60,0.85$を使います。

1. 作図。 床、点、法線、ランプを描き（1単位 = 2マス）、ランプの両端へのレイを引いて、法線との角度を測って下さい。
2. 正確な答え。 ランプは${\theta }_1$から${\theta }_2$を覆うので$E=L(\sin {\theta }_2-\sin {\theta }_1)$です。計算して下さい。
3. 一様な角度。 $\theta =\pi (u-\frac{1}{2})$でサンプリングします（pdf $1/\pi$）。各レイを描き、当たれば$L\cos \theta /(1/\pi )$、外れれば0です。4つの項を平均して下さい。
4. コサインサンプリング。 フラットランドのコサインのpdf $\frac{1}{2}\cos \theta$から逆関数法で$\theta =\arcsin (2u-1)$を導き、当たったサンプルは全て丁度$2L$を返す事を示し、4つの項を平均して下さい。
5. 光源サンプリング。 ランプ上の$s=\mathrm{\ell }u$を選び、$r$と$\cos \theta =\cos {\theta }_l=h/r$を求めます。フラットランドでは$\mathrm{d}\theta =\cos {\theta }_l\mathrm{d}s/r$なので$p_l(\theta )=r/(\mathrm{\ell }\cos {\theta }_l)$、各項は$L\cos \theta /p_l$です。レイを描き、項を平均して下さい。
6. 比較。 各手法の標本分散$\frac{1}{3}\sum (\text{項}-\text{平均}{)}^2$と、平均の誤差を求めて下さい。此処ではどれが一番良く、其のpdfの形でどう説明出来ますか？
7. MIS。 各光源サンプルについて$p_l$、$p_b=\frac{1}{2}\cos \theta$、$w_l=p_l/(p_l+p_b)$を求め、当たったコサインサンプルについては当たった点、$p_l$、$w_b$を求めて下さい。光源サンプル$k$とコサインサンプル$k$を組にすると、バランスヒューリスティックでは各項は$L\cos \theta /(p_l+p_b)$です。4つの組の推定値を平均して下さい。
8. 層化とQMC。 与えた数は$[0,1)$の4分の1毎に1つずつあり、層化された組です。ファン・デル・コルプット列の値$0.5,0.25,0.75,0.125,0.625,0.375,0.875,0.0625$で手順5をやり直し、4点の後と8点の後で平均して、誤差を比べて下さい。

自分で確かめる方法：正確な値は$E=\sin 45\mathrm{°}=0.7071$で、4サンプルでは平均は其の周りに散らばります（光源サンプリングは数%以内、一様とコサインは大きく外れる事もある）。
光源サンプリングの項はランプに沿って単調に減り、$L{h}^2\mathrm{\ell }/r_{max}^3=0.354$と$L\mathrm{\ell }/h=1$の間に入ります。
1つの方向の2つのMISの重みは足して1になり、此のランプの上では全て$w_l>\frac{1}{2}$です（$p_l=r^2/4\ge 1$、$p_b\le \frac{1}{2}$）。
終わった後の参考に、1サンプルあたりの正確な分散は1.519（一様）、0.914（コサイン）、0.0445（光源）です。