Pythonでパストレーサー

この記事は信州大学 kstmアドベントカレンダー21日目の記事です.

はじめに

皆様こんにちは.
誠に恐縮ながら21日目の記事を書かせていただきました. ~~よろですｗ~~

どんな記事が良いか考えたのですが,

光になりたい
レイトレ合宿[5]の作品見て感動した

という要求諸々があったので, 標題の内容になりました. 他にも以下のような内容が候補としてありました.

RasPiとカメラモジュールv2で品質(ISO9126, FURPS+等)を意識したデジカメを作る -> これはそのうち書く(多分)
論文執筆に便利な文献管理サービスを作った話 -> これもそのうち書く(多分)
VPNと過ごす素敵なゲーミング生活 -> これでも良いかなって思ったけど, なんか内容薄いかも？
実体験に基づく健康的なヘイトの稼ぎ方 -> 健康的などうかの評価が定量的に行えなかったので却下
◯◯から学ぶ言葉狩り -> 多分死ぬ
先輩を尊敬する文言集 -> 尊敬してます

などなど、書くと3日後くらいに僕の命が失くなりそうな内容のものは却下しました.

余談はさておき, 本題に入っていきます.

今回は, Pythonで物理ベースレンダリングに挑戦し, パストレーサーを作っていきます. 内容としては, hole氏が公開している[edupt][http://kagamin.net/hole/edupt/]をPythonで実装したものが主です. 正直, 中身の解説については解説PDFや文献[3]を見ていただいた方がまとまっていて分かりやすいと思いますが, 紹介や自分の勉強のために(アウトプット大事), 最低限で時間が許す限り書きたいと思います. (本当はGolangでやりたかった)

※数式について
既存の複数の文献ではベクトルxを\vec{x}と矢印を付けて表記しています. 私もこちらが分かりやすくて良いと思うんですが, 本記事では太字で\boldsymbol{x}と表記しています. 前者のほうが分かりやすいのでそのうち書き直すと思います.

レンダリング / Rendering

まず, レンダリングとは広義には以下のような定義です.[1]

レンダリング（rendering）とは、データ記述言語やデータ構造で記述された抽象的で高次の情報から、コンピュータのプログラムを用いて画像・映像・音声などを生成することをいう。

今回は写実的な綺麗な画像を生成することが目的ですので, 画像のレンダリングに着目します. 画像のレンダリングにおいて, 抽象的な高次の情報とは主にカメラの位置, 光の方向, オブジェクトの位置等がそれにあたります. (詳細は後述) 本記事では, それらの情報から物理モデルに基づくアルゴリズムでCGを生成することに着目します.

また, ゲームなどに応用される処理速度が重要なオンライン(リアルタイム)レンダリングではなく, 今回はそこまで目指さず, 処理時間がかかってもいいので綺麗な画像を作ることを目標とします.

物理ベースレンダリング / Physical Based Rendering

今回利用するレンダリング手法は, 光の物理モデルに基づく物理ベースレンダリングというものです. 英語でPysical Based Renderingというので, 頭文字を取ってPBRと略されます.

現実世界の光に関する物理現象を計算可能になるようにモデル化（数式で表現）して, それらを駆使することによって写実的な画像を得ることができます. 光に関する物理学は光学と呼ばれています.

光学(主に幾何光学)を基にしたレンダリング手法は, Ray(光線)の挙動のシミュレーションが主たる内容となっています. このシミュレーションを行うにあたって, 光のモデルに注目します.

以下に物理量について纏めた内容を載せますが, 自分も勉強中であることと, 参考文献[2,3,4,5,6]の方が綺麗に分かりやすく説明されているので, 正直勉強したいという方はそちらを見ることを推奨します.

放射束と放射照度 / Flux and Irradiance

単位時間あたりに面を通過するフォトン数（=エネルギー）を放射束Φ[W]と言い, 単位面積あたりの放射束を放射照度[W/m^2]と言います. 位置$x$, 面積$A[m^2]$に入る放射照度を$E(x)$とすると,
$$
E(x) = \frac{d\phi}{dA} \tag{1}
$$
で表されます.

立体角 / Steradian

普段使われる角度は平面角といい, 2次元上で定義されています. 一方立体角とはそれを3次元に拡張したものです. 立体角は球の半径を$r[m]$, 球上のある領域の面積を$A[m^2]$ とすると, 立体角$\Omega$は式(3)のように定義されます.
$$
\Omega = \frac{A}{r^2} \tag{2}
$$
証明等の解説は[こちらの動画][https://www.youtube.com/watch?v=aNoEzONgIYo]が一番分かりやすいと思います.(丸投げ)

放射輝度 / Radiance

微小立体角(単位立体角)あたりに到達する放射照度として定義されます.
$$
L(\boldsymbol{x},\boldsymbol{\omega}) = \frac{dE_{\omega}(\boldsymbol{x})}{d\omega}
$$
点$\boldsymbol{x}$において方向$\boldsymbol{\omega}$に通過するフォトン量(エネルギー)を放射輝度と捉える. Rayが物質を通過したり, 反射したりした時の物理量がこれにあたります. つまり, Rayが通過した時, ある点においてどのくらい放射輝度を発するかを計算することによってCGを生成することができます. 後に紹介するBRDF等を考慮することによって, より写実的なCGを生成することが可能になります.

レンダリング方程式

今回実装するレンダリングアルゴリズムでは式(2)の方程式を解きます.
$$
L_{o}(\boldsymbol{x},\boldsymbol{\omega}) = L_{e}(\boldsymbol{x},\boldsymbol{\omega})+\int_{\Omega_{x}}(f_{r}(\boldsymbol{x},\boldsymbol{\omega},\boldsymbol{\omega}')L_{i}(\boldsymbol{x},\boldsymbol{\omega})cos\theta) d \sigma(\boldsymbol{\omega}') \tag{2}
$$
$L_{O}(\boldsymbol{x},\boldsymbol{\omega'})$は, 点$\boldsymbol{x}$から方向$\boldsymbol{\omega}$への放射輝度を意味します. すなわち, Rayの反射や屈折などを経た後, どのくらいの輝度になるのかを示します. これを求めることによって最終的にカメラに入射する輝度を得ることができます.

$L_{e}(\boldsymbol{x},\boldsymbol{\omega})$は反射位置$\boldsymbol{x}$が光源上にある場合, 加算する放射輝度値です. 光源上にない場合は考慮しません.

$f_{r}(\boldsymbol{x},\boldsymbol{\omega},\boldsymbol{\omega'})$は$\boldsymbol{\omega'}$の方向から入射して$\boldsymbol{\omega}$方向へ進んでいく光がどの程度反射されるかを示す双方向反射率分布関数(Bidirectional Reflectance Distribution Function) と呼ばれます. 現実世界にはガラスや石など, 様々な物体が存在します. それらの物体の表面にはそれぞれ物体特有のザラつきや微小な穴などが存在します. これら全てを考慮するとあまりにも計算量が膨大になってしまいます. ですので, 物体がどの程度反射するかを(反射特性)をモデル化し, それに基づいて反射輝度を計算します. これによって計算量を抑えつつ物質特有の光り方をシミュレーションすることができます. また, 反射だけではなく透過や散乱についてのモデルも存在します(BTDF, BSDF)

$L_{i}(\boldsymbol{x},\boldsymbol{\omega'})$は位置$\boldsymbol{x}$に方向$\boldsymbol{\omega'}$から入射してくる放射輝度を意味します. 一つ前の反射によって得られた値がここに入ります.

$cos\theta$は, 平面に対して垂直に近いRayをより大きい放射輝度になるようにする項です. 平面に水平角度に入射してくるRayは考慮しないという解釈を私はしています.

これらの項を点$\boldsymbol{x}$を中心とした半球上の積分を計算することによって$L_{O}(\boldsymbol{x},\boldsymbol{\omega'})$を求めることが可能になります.
ただし, この積分方程式は$L_{i}(\boldsymbol{x},\boldsymbol{\omega'})$の項から分かる通り, 再帰的に解かなければなりません. このように複雑な積分方程式はモンテカルロ積分によって解きます.

モンテカルロ積分

まず1次元の定積分を考えます.　関数$f(x)$の$a$から$b$の区間の積分は
$$
\int_{a}^{b}f(x)dx
$$
と表されます. これをモンテカルロ積分で解こうとすると以下のように近似的に求めることが出来ます.
ただし, $N$はサンプル数, $x_{i}$は区間内で一様分布な確率変数です.
$$
\int_{a}^{b}f(x)dx \approx \frac{b- a }{N}\sum_{i=1}^{N}f(x_{i})
$$
時間の都合上話が飛躍してしまいますが, $x_{i}$を一様分布な確率変数ではなく, 任意の確率密度関数$p(x_{i})$に置換することが可能です. (ここの説明が曖昧なので多分書き直す)
$$
\frac{1}{N}\sum_{i=1}^{N}\frac{f(x_{i})}{p(x_{i})} + 誤差
$$
モンテカルロ積分では, 乱数を使って近似的に積分を求めます. これにより, 複雑な関数も近似的ですが積分を行うことが可能になります.

式(2)のレンダリング方程式をモンテカルロ積分によって解く場合, 式は以下のようになります.
$$
L_{O} = \frac{1}{N}\sum_{i=1}^{N}\frac{f_{r}(\boldsymbol{x},\boldsymbol{\omega},\boldsymbol{\omega}')L_{i}(\boldsymbol{x},\boldsymbol{\omega})cos\theta)}{p(\boldsymbol{\omega_{i}'})}
$$
$N回$乱数によって生成される$\boldsymbol{\omega_{i}'}$を元に, 放射輝度を計算し平均を取ることによって, 任意の点から出る放射輝度を近似的に求めることが可能になります. ($\boldsymbol{\omega_{i}'}$はある点を中心とした半球上を向いている)