まえがき
この記事は、42Tokyo Advent Calendar 2025の16日目の記事です。
私は今42Tokyoでレイトレーシングを使った課題に取り組んでいます。
今回用いたモンテカルロレイレーシング（パストレーシング）では、古典的なレイトレーシングに比べよりリアルな画像を作成できる一方、ノイズができやすくなってしまうという弱点があります。
今回はその回避方法の理論部分を学習しました。
図1: 純粋なモンテカルロレイトレーシング
図2: 今回の手法を用いたモンテカルロレイトレーシング
内容に誤り等あるかと思いますがご容赦ください。
ご指摘いただくと助かります。
参考資料
- 聖典
- 週末レイトレーシング 余生
- Ray Tracing: The Rest of Your Life(上の記事の最新版(英語))
- 三部作
- こちらもくわしい
- この記事を書くときに参考になった
モンテカルロ積分
確率密度関数
ある確率によって起こる事象に割り当てられる値を確率変数という。確率変数は離散型確率変数(サイコロの目など)と連続型確率変数(無作為に選んだ人の身長など)に分けられる。
ある確率変数の分布を縦軸に度数、横軸に階級をとってグラフにしたものをヒストグラムという度数を全項目数で割ればその階級が全体に占める割合(相対度数)を求められる。連続型確率変数についてこれを考えるのは無限に階級を分割することに等しい。ヒストグラムの階級を分割すると、その階級に含まれる度数が減る。もし無限に階級を分割すれば無限個の「度数0の階級」がうまれることになる。そこで連続型確率変数の場合は縦軸を各階級の度数ではなく全体に対する密度で表す。各階級の密度は以下の式で求められる。
階級の密度 =\frac{その階級にある項目の数}{全項目の数 ✕ 階級幅}
連続型確率変数においては確率変数 $X$ がある区間[a, b]の間に収まる確率 $P(a\leq X \leq b)$ を確率密度関数 $p(x)$ の積分によって以下のように表すことができる。
P(a\leq X \leq b) = \int_{a}^{b}p
(x)dx
したがって確率変数Xが微小区間[ $x, x+dx$ ]に含まれる値を取る確率は $p(x)dx$ である。
モンテカルロ積分とは
離散型確率変数Xの取りうる値を $x_1$, $x_2$ ,,, $x_n$ 、それぞれの確率を $p_1$, $p_2$ ,,, $p_n$ とすると、期待値Eは
E=\sum_{i = 1}^{n} x_ip_i
と表される。
連続型確率変数の場合、ある確率密度関数 $p(x)$ で分布する乱数 $x$ を用いると、確率が $p(x)dx$ であらわされることから、確率変数 $f(x)$ の期待値は
E = \int_{-∞}^{∞}f(x)p(x)dx
で表される。
これは、 $p(x)$ に従うN個の乱数 $x_1$ , $x_2$ ,,, $x_n$ を用いて、 $f(x_i)$ の平均を取ることで以下のように近似できる。
E = \int_{-∞}^{∞}f(x)p(x)dx ≈ \frac{1}{N}\sum_{i=1}^{N}f(x_i)
したがって、区間[a, b]におけるある関数f(x)の積分
I=\int_{a}^{b}f(x)dx
は
I=\int_{a}^{b}f(x)dx ≈ \frac{1}{N}\sum_{i=1}^{N}\frac{f(x_i)}{p(x_i)}
で近似できる。
大数の法則によれば，Nを十分に大きくすると，真値に収束すると考えることができる。
このように乱数を用いて積分結果の推定値を求める方法をモンテカルロ積分という。
累積分布関数
では、ある確率密度関数p(x)に沿った乱数を生成するにはどうすればよいだろうか。
それには以下のように定義される累積分布関数(cumulative distribution function, CDF) $P(x)$ という概念が必要になる。
P(x)=\int_{-∞}^{x}p(t)dt
CDFからPDFを求めるには導関数を取れば良い。
p(x) = \frac{d}{dx}P(x)
区間[0, 1]で一様分布する乱数 $u$ を用いて、とある確率密度関数 $p(x)$ に沿って乱数を生成する関数 $f(u)$ を考える。
$P(x)$ は $f(u)$ の値が $x$ よりも小さい確率を表す。
よってある定数 $a$ に対し
P_a = P(a)
が得られたとき、 $f(u)$ が $a$ より小さい確率は $P_a$ となる。ここで $u$ は0から1の間で一様に分布することを思い出す。 $f$ が単調増加すると仮定すると、 $P_a$ は $u$ の値が $0$ から $P_a$ の区間に入る確率と等しくなる。
したがって以下のように累積確率と出力 $f(u)$ の累積確率が対応すると考えられる。
f(P_a) = a
例えば $P(3) = 0.3$ つまり、 $f(u)$ の出力が $3$ より小さい確率が $0.3$ だったとすると、ランダムな入力 $u$ が $0.3$ を超えない確率は当然 $0.3$ であり $f(x)$ は入力が大きければ出力される値も大きくなるため $f(0.3) = 3$ と分かる。
したがって、以下のような関係が得られる。
f(P(x)) = x
すなわち $f$ は $P$ と逆の処理を行う。これを逆関数と呼び以下の様に表される。
f(x) = P^{-1}(x)
よって確率密度関数 $p(x)$ に従って分布す乱数を生成する関数 $f(x)$ は、その確率密度関数における累積分布関数 $P(x)$ の逆関数である。
重点サンプリング
いまいちど 先程示した下の積分を考える。
$x_i$ は区間[a, b]内で確率密度関数 $p(x)$ に従って分布しする。
I=\int_{a}^{b}f(x)dx ≈ \frac{1}{N}\sum_{i=1}^{N}\frac{f(x_i)}{p(x_i)} = \hat{I}_N
先程述べたようにモンテカルロ積分は、確率密度関数がpの乱数xを大量に生成し、それらに関する $\frac{f(x)}{p(x)}$ の平均を求める。この平均値が積分の結果の値に収束する。
期待値 $E[\hat{I}_N]$ は以下のように $f(x)$ の積分の結果と等しい。
E[\hat{I}_N] = \int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx = I
モンテカルロ法は高次元積分でも適用できる、計算が早いといった利点がある一方、実現値の収束が遅く、サンプリング数を大きくしないと真値に近づかないという欠点があるため、収束を早くする、つまり分散値を０に近づけることが重要になる。
分散 $V[\hat{I}_N]$ は
\begin{aligned}
V[\hat{I}_N] &= V\Biggr[\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}\Biggr] \\
&= \frac{1}{N^2}\sum_{i=1}^{N}V\Biggr[\frac{f(X_i)}{p(X_i)}\Biggr] \\
\end{aligned}
ここで
\begin{aligned}
V\Biggr[\frac{f(X)}{p(X)}\Biggr] &= E\Bigr[{\hat{I}_N}^2\Bigr] - E\Bigr[\hat{I}_N\Bigr]^2 \\
&= \int_{a}^{b}\Biggr(\frac{f(x)}{p(x)}\Biggr)^2p(x)dx - I^2 \\
&= \int_{a}^{b}\frac{f(x)^2}{p(x)}dx - I^2
\end{aligned}
$I$ は定数なので、 $\int_{a}^{b}\frac{f(x)^2}{p(x)}dx$ を最小化する $p$ を求めれば良い。ただしpは以下を満たす。
p(x)\geq0\ (a\leq x \leq b)\ ,\quad \int_{a}^{b}p(x)dx = 1
この値を小さくする方法は理論的に理解するには難解であるから、ここでは厳密な解析は割愛しもう少し感覚的に考えてみよう。
\int_{a}^{b}\frac{f(x)^2}{p(x)}dx
この値をできるだけ小さくするには、とにかく $p(x)$ を大きくすれば良いのではと思うかもしれないが $p$ は $\int_{a}^{b}p(x)dx = 1$ という条件があるため、 $p(x)$ の値は常に大きくはできない。どこかで $f(x)$ が $p(x)$ よりも著しく大きくなってしまう。ではどうするか。区間[a, b]で $p$ ができるだけ $f(x)$ と常に近い値をとれば良いのである。したがって $p$ が $f$ の近似として優れていればそれだけ収束も速くなる。 $p$ が $f$ と近いということはつまり、 $f(x)$ の値が大きい $x$ であるほど、サンプリングされる確率が高い確率分布を用いるということである。
球上の点を用いたモンテカルロ積分
先程は確率密度 $p(x)$ に基づいてサンプルする乱数 $x$ を用いて $f(x)$ を積分したが、我々がレイトレーシングで用いる乱数とは3次元空間の方向 $ω$ を意味している。方向 $ω$ は3次元空間の単位球面上の点として表すことができる。したがって方向 $ω$ による関数 $f(ω)$ の積分は球面上の面積 $S^2$ に関する積分として、3次元球面上の微小面積 $dσ$ を用いて次ように表される。
\int_{S^2}f(ω)dσ
ここでは例題として以下の積分を計算してみよう。
\int_{S^2}cos^2\theta dσ
方向 $ω$ は原点と単位球面上の点を結ぶ動径とz軸とのなす角 $\theta $ とxy平面への動径の射影とx軸のなす角 $\phi$ をもちいて、 $ω = (\theta , \phi)$ と表すことができる。これは半径が常に1の3次元極座標と言い換えてもよい。
単位球面の面積要素 $dσ$ には $dσ = sin\theta d\theta d\phi$ という等式が成り立つ。 (参照リンク)
したがって、
\begin{aligned}
\int_{S^2}cos^2\theta dσ &= \int_{0}^{2π}\int_{0}^{π}cos^2\theta sin\theta d\theta d\phi \\
\end{aligned}
ここで $u = cos\theta $ と置くと $du = -sin\theta d\theta $ である。
また、 $\theta = 0$ のとき $u = 1$ 、 $\theta = π$ のとき $u = -1$なので
\begin{aligned}
\int_{0}^{π}cos^2\theta sin\theta d\theta d\phi &= \int_{1}^{-1}u^2(-du) \\
&= \int_{-1}^{1}u^2du = \biggr[\frac{u^3}{3}\biggr]_{-1}^{1} = \frac{2}{3}
\end{aligned}
また
\int_{0}^{2π}d\phi = 2π
より
\begin{aligned}
\int_{0}^{2π}\int_{0}^{π}cos^2\theta sin\theta d\theta d\phi &= \Bigr(\int_{0}^{2π}d\phi\Bigr)\Bigr(\int_{0}^{π}cos^2\theta sin\theta d\theta d\phi\Bigr) \\
&= 2π⋅\frac{2}{3} = \frac{4π}{3}
\end{aligned}
となる。
これをモンテカルロ積分で近似的に求めるとどうなるだろうか。
確率密度 $p(ω)$ にしたがって方向 $ω$ をサンプリングし $cos^2\theta /p(ω)$ を計算する。ここでは方向 $ω$ を表す単位球面上の点が一様に選ばれるとすると $p(ω)$ は常に一定でかつ全体に渡る積分が1であるから、 $p(ω) = 1/(単位球面の面積) = 1/4π$ となる。また $cos^2\theta $ で $\theta $ はz軸との角度を表し、 $ω$ を極座標からxyz空間の直行座標へ変換すると $ω(\theta , \phi) = (sin\theta cos\phi, sin\theta sin\phi, cos\theta )$ のように表されるので、モンテカルロ積分の結果は次のように計算できる。
\begin{aligned}
\frac{1}{N}\sum_{i=1}^{N}\frac{cos^2\theta }{p(ω)} &= \frac{1}{N}\sum_{i=1}^{N}4πz^2
\end{aligned}
$z$ はランダムな単位ベクトルのz成分を表す。これをプログラムにしたのが以下である。
double pdf(direction p)
{
return (1 / (4.0 * M_PI));
}
#include <stdio.h>
int main(void)
{
int N = 147483647;
double sum = 0.0;
for (int i = 0; i < N; i++)
{
direction d = random_unit_direction(); // ランダムな方向を指す単位ベクトルを返す関数
double cos_squared = d.z * d.z;
sum += cos_squared / pdf(d);
}
printf("I = %f\n", sum / N);
}
// 出力例: I = 4.188611 (乱数を用いるため実行するたびに変化する。)
結果は $4π/3$ の近似になっている
光の散乱
今回行うレイトレーシングと積分にどんな関係があるのだろうか？
我々が行っているレイトレーシングとは、各ピクセルにレイを飛ばしそのピクセルの色を計算させる行為である。レイはカメラから放たれ、オブジェクトとの衝突点で散乱しまた次のレイを飛ばす。したがってレイトレーシングにおいて色の計算ではレイを引数として色を返す関数を作成し、衝突点を中心とする半球上の点で表されるあらゆる方向についてこの関数を積分する。そしてそれを拡散したレイについても再帰的に繰り返すことでシュミレーションを行うのである。
ランバート反射
ここで光度と輝度の違いについて理解しおこう。
光度とはある方向へ放たれる光の強さを表す量である。一方、輝度とは単位面積あたりの光度を表す量である。
ある反射面の微小領域 $dA$ がその法線 $N$ から角度 $\theta $ 分斜めの向きで観測されたとき、 $dA$ が $N$ の向きに放つ光度 $dI$と観測者が受け取る輝度 $L$ の関係は以下の式で表される。
L = \frac{dI}{dA\cdot cos\theta }
どの角度から見ても均一な明るさで見える理想的な拡散反射をランバート反射という。ランバート反射をする表面はランバート面と呼ばれる。
ランバート面では輝度が常に一定になるため、上の式より光度 $I$ は $cos\theta $ に比例する。
アルベド
物体に入射した光のうちすべてが反射されるわけではなく一部は物体に吸収される。プログラムでは物質ごとの光を反射する割合 $A$ で表しコレをアルベドと呼ぶ。反射率は色によって異なるためRGBの3値で表される。光が吸収される確率は $1 - A$ で表すことができる。
散乱
光源から放たれた光は物体との衝突点で、物質の反射率の影響を受けながら散乱しカメラのレンズや私達の水晶体に到達することで色として認識される。レイトレーシングでは実際の物理現象とは逆にカメラからレイを放ち衝突点に入射してきた無数の光をさらにランダムにサンプリングした大量のレイを飛ばして逆向きに追跡していくことで色を計算する。
立体角
球面上のある部分の面積に対し球の中心からどの程度の広がりを持つかを表す量を立体角という。立体角 $Ω$ は半径 $r$ の球面上のある部分の面積を $S$ として時以下の式で求められる。
Ω = \frac{S}{r^2}
単位球においては、立体角の大きさは球面上の面積と等しい。
散乱PDF
光が衝突点から散乱するときの各方向の光度を散乱レイの分布で表すと、散乱する方向 $ω_o$ を引数にとるPDFとして表される。これを以降は散乱PDFと呼び、 $pScatter(ω_o)$ と表す。散乱PDFをある立体角で積分するとその立体角の内側にレイが散乱する確率を求めることができる。散乱PDFは入射方向によっても変化する場合があるため次のように表す: $pScatter(ω_o, ω_i)$ 。
物体の表面で散乱する光の色 $Color_o$ は入射してくる光の色 $Color_i$ とアルベド $A$ 、散乱PDF $pScatter(ω_o, ω_i)$ の積をすべての入射方向 $Ω_i$ の範囲で積分することで次のように表せる。
Color_o = \int_{Ω_i}Color_i\cdot A\cdot pScatter(ω_o, ω_i)dω_i
入射光の色 $Color_i$ もまたその光が反射してきた物体表面で同様の計算をして求められるためこれは再帰的なアルゴリズムである。
ランバート面の散乱PDF
先程の積分を、モンテカルロ積分を使って推定していきたい。一回の試行では入射方向 $ω_i$ に対してその入射方向のPDF $p(ω_i)$ を用いて、以下の式を計算することになる。
Color_ o = \frac{Color_i\cdot A\cdot pScatter(ω_o, ω_i)}{p(ω_i)}
ここでランバート反射の $pScatter(ω_o, ω_i)$ を考える。前述したように物体表面での法線と散乱レイの角度を $\theta $ とするとランバート面での光度(ここでは $pScatter$ )は $cos\theta $ に比例する。つまり $pScatter(ω_o, ω_i) = C\cdot cos\theta $ と表せる。 $cos\theta <0$ つまり散乱レイが表面より内側に散乱することはないので $pScatter(π/2<0\leq π) = 0$ とする。 またここでPDFの積分は1になることを思い出す。半球 $Ω_i$ に関して $pScatter$ の積分し
\begin{aligned}
\int_{Ω_i}pScatter(ω_o, ω_i)dσ &= \int_{Ω_i}C\cdot cos\theta dσ \\
&= \int_{0}^{2π}\int_{0}^{\frac{π}{2}}cos\theta sin\theta d\theta d\phi\cdot C \\
&= 2π\cdot \frac{1}{2}\cdot C \\
1 &= π\cdot C \\
C &= \frac{1}{π}
\end{aligned}
となることから
pScatter(ω_o, ω_i) = \frac{cos\theta }{π}
である。
散乱PDFと等しいPDFを使ってサンプルを行うとすれば
p(ω_i) = pScatter(ω_o, ω_i) = \frac{cos\theta }{π}
より今節先頭の等式の分母と分子が打ち消し合い
\begin{aligned}
Color_ o &= A\cdot Color_i
\end{aligned}
が得られる。 これは、
t_color ray_color( /* ... */ )
{
// ...
return (albdo * ray_color( /* ... */ ));
}
に等しい。
重点サンプリング付きマテリアル
いまいちど以下の式みてみよう。
Color_ o = \frac{Color_i\cdot A\cdot pScatter(ω_o, ω_i)}{p(ω_i)}
先程はランバート反射における $pScatter(ω_o, ω_i)$ と等しい確率密度関数 $p(ω_i)$ を用いた。したがってサンプリングされる方向 $ω_i$ の分布は $cos\theta $ に比例する。
モンテカルロ積分においては、被積分関数の値が大きくなる乱数を多くサンプリングすることで分散が小さくなり収束を早くすることができるのであった。レイトレーシングにおいては収束が早くなることはすなわち、同数のサンプル数でよりノイズが少なくなることを意味する。
今回の被積分関数とは以下の部分で
Color_i\cdot A\cdot pScatter(ω_o, ω_i)
この値が大きいとはつまり、より明るいということである。
これから画像のノイズを減らすためにレイを光源に向かって放つ確率が高くなるようにすることを考える。光源へ多くレイを放つときのPDFを $pLight(ω_i)$ として、 $pScatter(ω_i)$ が大きい部分へレイを放つときのPDFをこれからは $pSurface(ω_i)$ と呼ぶことにする。
PDFには複数のPDFの線形和を取って確率密度を混ぜることで別の PDF を作れるという便利な性質がある。例えば
p(ω_i) = \frac{1}{2}pSurface(ω_o, ω_i) + \frac{1}{2}pLight(ω_o, ω_i)
以上の $p(ω_i)$ はPDFの性質を満たしている。
先程述べたように我々の目的は以下の積
Color_i\cdot pScatter(ω_o, ω_i)
が大きい値となるようなPDF $p(ω_i)$ を見つけることである。( $A$ は定数なので除いて良い)
ランバート面の場合、 $pScatter(ω_i)$ は入射角に依らないため、 $Color_i$ が大きくなる、つまりより明るい入射方向 $ω_i$ を多くサンプリングすれば良い。
鏡面の場合、 $pScatter(ω_i)$ が特定の方向の周りで非常に大きくなるため $pScatter(ω_i)$ が大きくなる方向をより多くサンプリングする必要がある。このため多くのレンダラでは鏡を例外的なケースとして、以下のように扱っていることが多い。
Color_ o = \frac{Color_i\cdot A\cdot pScatter(ω_o, ω_i)}{p(ω_i)} ≈ \frac{pScatter(ω_i)}{p(ω_i)}
ランダムな方向の生成
何らかの分布に沿った方向 (単位球上の点) をランダムに生成する方法を見つけよう。ここでは逆関数法というやり方を用いてみる。
まずは簡単のため z 軸を法線方向として、 $\theta $ を法線から測った角度とする。軸を回転させる方法は次節で説明する。考える分布はz軸周りの回転に関して対称な分布だけとする。このときPDFは $\theta $ の関数で表せる。
方向を引数に取るPDF $p(ω) = f(\theta )$ が与えられたとして、 $\theta $ 及び $\phi$ に関する一次元のPDFを $a(\phi)$ , $b(\theta )$ とすると、 $\theta $ と $\phi$ は独立した確率変数なので $p(ω) = f(\theta ) = a(\phi)\cdot b(\theta ) $ と表せる。微小面積 $dσ$ がサンプルされる確率は
p(ω)dσ = b(\theta )a(\phi)d\theta d\phi
また $dσ = sin\theta d\theta d\phi$ より
p(ω)dσ = f(\theta )sin(\theta )d\theta d\phi
$a(\phi)$ は区間[ $0$ , $2π$ ]で一様なので $a(\phi) = 1/2π$ 。
よって
\begin{aligned}
b(\theta )a(\phi)d\theta d\phi &= f(\theta )sin(\theta )d\theta d\phi \\
b(\theta )a(\phi) &= f(\theta )sin(\theta ) \\
b(\theta )\cdot \frac{1}{2π} &= f(\theta )sin(\theta ) \\
b(\theta ) &= 2πf(\theta )sin\theta
\end{aligned}
したがって
\begin{aligned}
a(\phi) &= \frac{1}{2π} \\
b(\theta ) &= 2πf(\theta )sin\theta
\end{aligned}
が成り立つ。
区間[0, 1]の一様乱数 $u_1$ を用いて、 $a(\phi)$ にしたがうランダムな角度 $\phi$ を生成するには、 $a(\phi)$ の累積密度関数の逆関数を求める。これには関数の入力 $u_1$ と出力 $\phi$ を反転し累積密度関数 $P(x)$ に与えた式を考えれば良い。
\begin{aligned}
u_1 &= P(\phi) \\
&= \int_{0}^{\phi}a(t)dt = \int_{0}^{\phi}\frac{1}{2π}dt \\
&= \frac{\phi}{2π}
\end{aligned}
これを $\phi$ について解くと
\phi = 2π\cdot u_1
を得る。これで $u_1$ を用いて $a(\phi)$ に沿ったランダムな角度 $\phi$ を生成できる。
同様に区間[0, 1]の一様乱数 $u_2$ に対して $b(\theta )$ の累積密度関数の逆関数を求めると
\begin{aligned}
u_2 &= \int_{0}^{\theta }b(t)dt \\
&= \int_{0}^{\theta }2πf(t)sin\ t\ dt
\end{aligned}
ここで $f$ はいまランダムな方向のPDFを表している。仮にここではすべての球面上の点に一様な分布を考えよう。単位球の表面積は $4π$ だから一様分布では $p(ω) = f(\theta ) = 1/4π$ となるので
\begin{aligned}
u_2 &= \int_{0}^{\theta }2πf(t)sin\ t\ dt \\
&= \int_{0}^{\theta }2π\cdot \frac{1}{4π} sin\ t\ dt \\
&= \int_{0}^{\theta }\frac{1}{2} sin\ t\ dt \\
&= \frac{-cos\theta }{2} - \frac{-cos(0)}{2} \\
&= \frac{1 - cos\theta }{2}
\end{aligned}
であり $cos\theta $ について解くと
cos\theta = 1 - 2u_2
を得る。 $\theta$ より先に $cos\theta$ が求まる場合が多いから、この式を $\theta$ について解く必要はない。これで $b(\theta )$ に沿ったランダムな $cos\theta $ を生成できる。
方向 ( $\theta$ , $\phi$ ) に向かう単位ベクトルを生成するには、次の式を使ってデカルト座標系に変換する必要がある
\begin{aligned}
x &= cos\phi\cdot sin\theta \\
y &= sin\phi\cdot sin\theta \\
z &= cos\theta
\end{aligned}
これと恒等式 $cos^2x+sin^2x=1$ を使えば、二つの一様乱数 $u_1$ , $u_2$ を使った表現が求まる。
\begin{aligned}
x &= cos(2π\cdot u_1)\sqrt{1-(1-2u_2)^2} \\
y &= sin(2π\cdot u_1)\sqrt{1-(1-2u_2)^2} \\
z &= 1 - 2u_2
\end{aligned}
$(1 - 2u_2)^2 = 1 - 4u_2 + 4(u_2)^2$ より
\begin{aligned}
x &= cos(2πu_1)\cdot 2\sqrt{u_2(1-u_2)} \\
y &= sin(2πu_1)\cdot 2\sqrt{u_2(1-u_2)} \\
z &= 1 - 2u_2
\end{aligned}
となる。これを使って点を生成すると以下のようになる。
#include <math.h>
int main(void)
{
for (int i = 0; i < 200; i++)
{
double u1 = random_double(0.0, 1.0);
double u2 = random_double(0.0, 1.0);
double x = cos(2*M_PI*u1)*2*sqrt(u2*(1-u2));
double y = sin(2*M_PI*u1)*2*sqrt(u2*(1-u2));
double z = 1 - 2*u2;
printf("%f %f %f\n", x, y, z);
}
return (0);
}
正規直行基底
前節では z 軸を基準にランダムな方向を生成する方法を導出した。次は物体表面の法線を基準に同じことを行いたい。
相対座標
互いに直交する三つの単位ベクトルの集合を正規直交基底 (OrthoNormal Basis) と呼ぶ。例えばデカルト座標系の x, y, z 軸方向の単位ベクトルからなる集合は ONB である。
現実世界における物体の位置と傾きを決めるには、基準となる ONB がどこかに必要となる。写真とはカメラから見たシーンの相対的な位置と傾きを写したものであり、カメラとシーンが同じ座標系で表されている限りどんな座標系でも問題はない。
正規直交基底の構築
法線 $\vec{n}$ を基準としたONBを構築する。長さが $0$ より大きく $\vec{n}$ と平行でない適当なベクトルを $\vec{a}$ とすれば、 $\vec{n}$ と垂直で互いに直交する二つのベクトル $\vec{s}$ , $\vec{t}$ は次の式から得られる。
\begin{aligned}
\vec{t} &= unit\_ vector(\vec{a} × \vec{n}) \\
\vec{s} &= \vec{t} × \vec{n}
\end{aligned}
よく考えずに適当な $\vec{a}$ を選ぶと、 $\vec{n}$ に平行な $\vec{a}$ を選んでしまう可能性がある。
これに対処するためによく使われるのが、if文 を使って $\vec{n}$ がxyz軸のいずれか平行かどうかを判定し、もし平行でなければその軸を $\vec{a}$ として採用するという方法である
n = normalize(n);
if absolute(n.x > 0.9)
a = (0, 1, 0)
else
a = (1, 0, 0)
ONBが計算できたら、次に z 軸を基準にしたランダムな方向ベクトル $random(x,y,z)$ を生成する。すると法線 $\vec{n}$ を基準にしたランダムな方向ベクトルは
x\vec{s} + y\vec{t}+ z\vec{n}
と表せる。
直接的なライトのサンプリング
ほぼ一様に方向をサンプルする方法の問題は、ライトが存在する方向と他の重要でない方向が同程度にしかサンプルされないことである。ここではライトのある方向へ多くレイを放つようにする。
ライトに向かうランダムな方向を選ぶのは非常に簡単で、ライト上の点をランダムに選び、その点に向かう方向を選べばよい。ただモンテカルロ積分を行うには、その方向の PDF p(direction) を知る必要がある。PDF はいくつだろうか？
ライトのPDF
まずは平面の光源を考えよう。面積 $A$ のライト上の点を一様ランダムに選択すると、ライト上の任意の点のPDFは $1/A$ となる。方向を定義する単位球上におけるライトの見かけの面積(立体角)を求めよう。ライト上に小面積 $dA$ を取り、その内部の点を適当に $q$ とする。このとき一様ランダムに選択したライト上の点がこの小面積内にある確率は $dA/A$ である。また求めようとしている PDF $p(ω)$ において立体角で表される単位球上の小面積 $dσ$ をサンプルする確率は $p(ω)⋅dσ$ と表せる。さらに、単位点の中心点 $p$ とライト表面の点 $q$ の距離を $distance(p, q)$ とすると、立体角の定義より $dσ$ と $dA$ の間には次の関係がある。
dσ = \frac{dA\cdot cos\alpha}{distance^2(p, q)}
ここで $\alpha$ は $ω$ とライトの法線がなす角度を表す。 $dσ$ と $dA$ から点が選ばれる確率は同じだから
p(ω)\cdot dσ = p(ω)\cdot \frac{dA\cdot cos\alpha}{distance^2(p, q)} = \frac{dA}{A}
つまり $0\leq\alpha< π/2$ において
p(ω) = \frac{distance^2(p, q)}{cos\ \alpha\cdot A}
が成り立つ。( $π/2\leq\alpha< π においてp(ω) = 0$ )
これはライトに向かう立体角の中から1つの方向が選ばれるとしたときの確率密度を表しているためライトの見かけの面積が小さければ確率分布は集中するため $p(ω)$ は大きくなる。
球オブジェクトのサンプリング
球の外側から球の表面を立体角に関して一様にサンプルするとき、球に接する円錐をサンプルすると考えても同じことになる 。始点 $p$ から見える球の表面上の点を一様にサンプリングすることは、円錐の軸と方向のなす角 $\theta $ と、円錐の軸を基準とした一様な方位角 $\phi$ をサンプリングすることに等しい。
半径 $R$ の球の中心 $c$ と $p$ を結ぶ直線(円錐の軸)と、$p$ を通り球と接する直線のなす角を $\theta _{max}$ とすると $\theta $ は $0\leq \theta \leq \theta _{max}$ の範囲で一様に選ぶことになる。この一様なサンプリングのPDFは定数 $C$ と表される。区間[0, 1]を用いて乱数を生成する関数を先ほどと同様に逆関数法で求める。
\begin{aligned}
u_1 &= \int_{0}^{\theta }2π\cdot C\cdot sin\ t\ dt \\
&= 2π\cdot C\cdot(1 - cos\theta ) \\
cos\theta &= 1 - \frac{u_1}{2π\cdot C}
\end{aligned}
を得る。 $u_1 = 1$ で $\theta = \theta _{max}$ より
\begin{aligned}
cos\theta _{max} &= 1 - \frac{1}{2π\cdot C} \\
\frac{1}{C} &= 2π(1 - cos\theta _{max})\\
\end{aligned}
よって
cos\theta = 1 + u_1(cos\theta _{max} - 1)
を得る。 $\phi$ は $0 \leq \phi < 2π$ で一様なので
\phi=2π⋅u_2
で求められ、PDFは $1/2π$ である。
よって球に向かう方向 $ω$ は
\begin{aligned}
cos\theta &= 1 + u_1(cos\theta _{max} - 1) \\
\phi &= 2π⋅u_2
\end{aligned}
を満たす $ω(\theta , \phi)$ で生成することができ、これを直行座標系に変換すると
\begin{aligned}
x &= cos\phi\cdot sin\theta = cos(2π⋅u_2)\cdot \sqrt{1 - z^2}\\
y &= sin\phi\cdot sin\theta = sin(2π⋅u_2)\cdot \sqrt{1 - z^2} \\
z &= cos\theta \hspace{10mm} = 1 + u_1(cos\theta _{max} - 1)
\end{aligned}
となる。
続いて $\theta_{max}$ を考える。
$\sin\theta_{max} = R / distance(c, p)$ より
\cos\theta_{max} = \sqrt{1 - \frac{R^2}{distance^2(c, p)}}
である。そして方向のPDFは $1/(球を覗く立体角)$ と表せるので
\begin{aligned}
球の立体角 &= \int_{0}^{2\pi}\int_{0}^{\theta_{max}}\sin\theta d\theta d\phi \\
&= 2\pi(1 - \cos\theta_{max})
\end{aligned}
より、PDFは $1/2\pi(1 - \cos\theta_{max})$ と求まる。
混合密度
これで表面で拡散するときの $\cos\theta$ に比例するPDFとライトを直接サンプルするPDFが手に入った。この２つを組み合わせたPDFを作りたい。
確率論では複数の密度関数を混ぜて混合密度を作るテクニックがよく使われる。例えば二つの密度関数の平均を取れば混合密度となる
p_{mixture}(ω) = \frac{1}{2}p_{surface}(ω) + \frac{1}{2}p_{light}(ω)
レイトレーサーのアーキテクチャについて
・MISを利用するとシャドウレイは使えなくなる。
・鏡面反射と拡散反射どちらも行うような表面の場合は $p_{surface}$ を２通り用意し確率的に選択すれば良い。
・また完全な鏡面反射を実装する際には例外的に処理を分ける必要がある。（粗さが0でない鏡面反射は実装できる）
・本格的な物理ベースレンダラではRGBではなく光のスペクトルをもとにしたパラメーターで色を扱う。