物理ベースレンダラ

edupt解説
2014/06/13 Ver.1.03
hole (@h013)
http://kagamin.net/hole/edupt/index.htm

eduptとはシンプルでコンパクトな物理ベースレンダラ

eduptとはsmallpt をベースに修正・拡張

C++で記述日本語によるコメント付き

教育的

smallpt (http://www.kevinbeason.com/smallpt/) はKevin Beason氏に
よってC++で書かれた99行のパストレーシングによるレンダラ

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

このスライドは

eduptのコード解説＋
基礎的な物理ベースCGの解説＋
基礎的なパストレーシングの解説

目次
1.
2.

物理ベースコンピュータグラフィックス
eduptコード解説（ユーティリティ編）
1.
2.
3.

3.

4.

5.

各種定数
画像出力
乱数

eduptコード解説（幾何編）
1.
2.
3.

ベクトル
球
シーンデータ

1.
2.
3.
4.
5.

光の物理量
レンダリング方程式
モンテカルロ積分
パストレーシング
ロシアンルーレット

パストレーシング

eduptコード解説（レンダリング編）
1.
2.
3.
4.
5.

カメラ設定
画像生成
radiance()
完全拡散面
完全鏡面・ガラス面

1.物理ベースコンピュータグラフィックス

物理ベースコンピュータグラフィックス まえおきeduptはパストレーシングと呼ばれるアルゴリズムによる物理ベースレンダラですが、そもそも物理ベースCGとは何なのか？何が嬉しいのか？何ができるのか？について

物理ベースコンピュータグラフィックス

何らかの物理モデルに基づいて計算されたコンピュータグラフィックス

物理モデル現実世界における各種物理現象を数式で表現したもの• 物理モデルによって各種物理現象が近似的に計算可能になる。

𝑚𝑎 = 𝐹ニュートン力学も物理モデルの一種

物理ベースじゃないかもしれないコン
ピュータグラフィックス
 何らかの物理モデルに基づいていても、十分にそのモデルに従っていない計算によって得られたCG• 光の相互反射を考慮していないようなCG

 そもそも何の物理モデルに基づいていない計算によって得られたCG• エネルギー保存則に従っておらず、光がどんどん増幅していくようなCG物理ベースの嬉しいところ
 現実の近似としての物理モデルを使う＝現実らしいCGが得られる！（＝フォトリアル）
 各種パラメータ（材質や光源の設定等）について、現実世界の値をそのまま使える！
 既に存在する物理モデルを基礎においているため、理論的な安心感がある！
 などなどコーネルボックス

光学どんな物理モデルを使う？
CGにおいて考慮すべき物理モデルは光に関する物理学、光学についてのモデル
様々な光学モデル• 幾何光学• 波動光学• 量子光学• etc

現代のCGでは、基本的に幾何光学をモデルとして用いる。これは、目に映る風景や写真によって描き出される像が幾何光学によって十分再現できるからである。
※もちろん、幾何光学では再現できない物理現象もあり、そういった現象を再現したい場合は別のモデルを採用する。

幾何光学
 光の伝播を光線として記述• 光線＝Ray• いわゆるレイトレーシングのレイとはこの光線のこと。

CG的には何ができるのか？
 幾何光学で表現できることは何でもできる• グローバルイルミネーション• 被写界深度• モーションブラー
• etc
eduptでは特に光の相互反射による影響を考慮した
グローバルイルミネーションの実現を重視している。幾何光学でグローバルイルミネーションを実現するための物理モデルはこのスライドの後々取り扱う。

2.eduptコード解説（ユーティリティ編）

2.eduptコード解説（ユーティリティ編）

 まえおきeduptにおける、各種ユーティリティ関数についての解説。

eduptファイル概要

main.cpp
edupt::render()を呼び出すだけ。

constant.hintersection.hmaterial.h

各種定数。
交差判定の結果を格納する構造体IntersectionとHitpoint。
Color型と物体材質について記述するReflectionTypeと屈折率の設定。



ppm.h



radiance.hrandom.h



ray.h



render.hscene.h

sphere.hvec.h

レンダリング結果をppm画像として書き出すための関数save_ppm_file()。
ある方向からの放射輝度を得る関数radiance()。パストレーシングのメイン処理。
乱数生成クラスXorShift。
一つ一つの光線、レイを表現する構造体Ray。
レンダリング画像サイズやサンプル数を受け取り、radiance()を使って各ピクセルの具体的な値を決定する関数
render()。
レンダリングするシーンのデータspheres[]とそのシーンに対する交差判定を行う関数intersect_scene()。

基本的な形状としての球を表現するSphere構造体。
ベクトルを表現する構造体Vec。

eduptファイル概要

main.cpp
edupt::render()を呼び出すだけ。

constant.hintersection.hmaterial.h

各種定数。
交差判定の結果を格納する構造体IntersectionとHitpoint。
Color型と物体材質について記述するReflectionTypeと屈折率の設定。



ppm.h



radiance.hrandom.h



ray.h



render.hscene.h

sphere.hvec.h

レンダリング結果をppm画像として書き出すための関数save_ppm_file()。
ある方向からの放射輝度を得る関数radiance()。パストレーシングのメイン処理。
乱数生成クラスXorShift。
一つ一つの光線、レイを表現する構造体Ray。
レンダリング画像サイズやサンプル数を受け取り、radiance()を使って各ピクセルの具体的な値を決定する関数
render()。
レンダリングするシーンのデータspheres[]とそのシーンに対する交差判定を行う関数intersect_scene()。

基本的な形状としての球を表現するSphere構造体。
ベクトルを表現する構造体Vec。

2.1 各種定数

constant.h 各種定数を記録したヘッダ– kPI• 円周率π。

– kINF• 非常に大きい定数。

– kEPS• 0.0に非常に近い定数。交差判定等で使う。（球の交差判
定のところで解説）2.2 画像出力

ppm.hsave_ppm_file() レンダリングした結果を保存するための関数– Color型の配列を受け取る。
– Color型はRGBデータをそれぞれx,y,zメンバ変数で保存している。

save_ppm_file() eduptはppm形式で保存– http://netpbm.sourceforge.net/doc/ppm.html– ppmはテキスト形式で画像を記録できるフォーマット。
– ヘッダに”P3 (改行) 幅 高さ (改行) 255”と記録することで、指定したサイズのRGB画像（RGBの階調は256段階）を記録する、という指定になる。
– ヘッダに続けて”R値 G値 B値“を画素数の分だけ並べる。

clamp()とto_int() レンダリング結果の配列（Color型）の各値はto_int関数で0から255の256階調の値に変換されて記録される
1. 各値xはclamp()によって0から1の範囲に詰められる
2.

1
乗する（ガンマ補正）
2.2

3. 255倍する（0から255の範囲になる）
4. 0.5を足してintにキャストする（小数第一位で四捨五入）

ガンマディスプレイのガンマ値– 一般的なディスプレイは入力値と出力値の間に非線形な関係が成り立つ。
プログラム側の輝度値（入力）

ディスプレイによって出力される
物理的な輝度値

1.0
0.7
0.5
0.3
0.0

1.0
0.46
0.22
0.07
0.0

出力= 入力

𝛾

ガンマ値

ガンマ
 γの値
– 基本的にディスプレイ、OSに依存– sRGB、AdobeRGBのような広く業務にも使われてるカ
ラースペースではγ=2.2が採用されている。
1.0

γ=1.0

出力 0.5

γ=2.2

0.0
0.0

0.5
入力

1.0

ガンマ補正 物理ベースレンダリングの結果はリニアスペースなのでディスプレイの出力輝度値を画像の値に比例させたい– ガンマ補正– γ値の逆数のべき乗をディスプレイに対する入力にする。

– eduptではγ=2.2と仮定して画像の保存の際値を1乗している。
2.2

プログラム側の輝度値
（入力）

補正値

ディスプレイによって出力される
物理的な輝度値

1.0
0.7
0.5
0.3
0.0

1.0
0.85
0.73
0.58
0.0

1.0
0.7
0.5
0.3
0.0

ガンマ補正

補正値= 入力1𝛾𝛾

出力= 補正値

ガンマ補正ガンマ補正前ガンマ補正後γ=2.2と仮定

2.3 乱数

random.h乱数
 パストレーシングは確率的なアルゴリズム– さまざまな局面で乱数が必要になる。
– 質の高い乱数が求められる。

 擬似乱数アルゴリズム– 線形合同法– メルセンヌツイスター– Xorshift– Random123– etc

eduptでは実装の単純さと速度、質などのバランスからXorshiftを乱数生成機として使用している

パストレーシングがモンテカルロ積分をその基礎においているため、サンプルを得るために乱数が必要になる。準モンテカルロ法の様な、より進んだアルゴリズ
ムを使うなら、乱数の変わりに低食い違い量列というものを使うこともある。

Xorshift eduptでの実装– WikipediaのXorshiftの項目• http://ja.wikipedia.org/wiki/Xorshift– Xorshiftのseed初期化• http://meme.biology.tohoku.ac.jp/klabowiki/index.php?%B7%D7%BB%BB%B5%A1%2FC%2B%2B– を参考にしました。

– next01()によって[0.0, 1.0]の範囲で乱数を得る。

 Random型にtypedefして他の場所で使用

3.eduptコード解説（幾何編）

3.eduptコード解説（幾何編）

 まえおきeduptにおける、幾何学関係のデータ構造についての解説。
ベクトル構造体や、基本的な形状である球とレイの交差判定など。

eduptファイル概要

main.cpp
edupt::render()を呼び出すだけ。

constant.hintersection.hmaterial.h

各種定数。
交差判定の結果を格納する構造体IntersectionとHitpoint。
Color型と物体材質について記述するReflectionTypeと屈折率の設定。



ppm.h



radiance.hrandom.h



ray.h



render.hscene.h

sphere.hvec.h

レンダリング結果をppm画像として書き出すための関数save_ppm_file()。
ある方向からの放射輝度を得る関数radiance()。パストレーシングのメイン処理。
乱数生成クラスXorShift。
一つ一つの光線、レイを表現する構造体Ray。
レンダリング画像サイズやサンプル数を受け取り、radiance()を使って各ピクセルの具体的な値を決定する関数
render()。
レンダリングするシーンのデータspheres[]とそのシーンに対する交差判定を行う関数intersect_scene()。

基本的な形状としての球を表現するSphere構造体。
ベクトルを表現する構造体Vec。

3.1 ベクトルvec.hVec構造体
三組の実数による構造体– 線形代数におけるベクトルとほぼ同じ。
– eduptでは点、ベクトル、色等を表現するために使われる。
ｙｚｘ
三次元空間におけるベクトル

Vec構造体 ベクトル同士の加算
 ベクトル同士の減算 ベクトルとスカラ量の乗算 ベクトルとスカラ量の除算
 二乗長さ、長さを求める関数length_squared()とlength()– 各要素の二乗の和の平方根がベクトルの長さ（三平方の定理）

normalize()ベクトルの正規化関数– ベクトルの各要素を自身の長さで割る操作。
– 正規化後のベクトルは長さが１になる。
𝑣=
𝑣

𝑣
𝑣𝑥2 + 𝑣𝑦2 + 𝑣𝑧2

multiply()二つのベクトルの要素ごとの積dot()
二つのベクトルの内積を計算する
内積𝑢∙𝑣 = 𝑢𝑣 cos 𝜃𝑢 ∙ 𝑣 = 𝑢𝑥 𝑣𝑥 + 𝑢𝑦 𝑣𝑦 + 𝑢𝑧 𝑣𝑧

𝜃𝑢𝑣

特に長さ１のベクトル同士の内積は二つのベクトルの成す角のcosとなる。
二つのベクトルの成す角を得たいときは内積を取るのが定石。

cross()二つのベクトルの外積を計算する
外積

𝑢 × 𝑣 = (𝑢𝑦 𝑣𝑧 − 𝑢𝑧 𝑣𝑦 , 𝑢𝑧 𝑣𝑥 − 𝑢𝑥 𝑣𝑧 , 𝑢𝑥 𝑣𝑦 − 𝑢𝑦 𝑣𝑥 )外積によって、二つのベクトルのどちらとも直交するベクトルが得られる。
その向きは右ネジの進む方向になる。

外積の長さ

𝑢×𝑣 = 𝑢𝑣 sin 𝜃𝑢𝜃

外積の長さは二つのベクトルの長さと成す角のsinとの積になる。
また、二つのベクトルによって作られる平行四辺形の面積に等しい。

𝑢×𝑣𝑣

ray.h𝑑𝑖𝑟 レイ（光線）を表現する構造体
𝑜𝑟𝑔– 始点orgと向きdirを持つ。
– 始点は位置ベクトルとしてVecで表現する。
– eduptでは基本的にdirの長さは1になるようにしている。

• 何か所かこの仮定を用いている処理がある。
• もちろんこの仮定をなくしてもよい。（ただしこの仮定を用いている処理は修正しないといけない）

3.2 球sphere.hSphere構造体
 eduptにおける唯一の幾何形状– 中心座標 positionと半径 radius– 発光色 emission– 反射率 color– 材質（表面における反射の種類）reflection_type

radiuspositionmaterial.h Vec構造体を使ってColorを定義
 材質ReflectionTypeと屈折率– 完全拡散面、完全鏡面、ガラス面の三種類– ガラスに対する屈折率がkIor– 詳しくはレンダリング編

球との交差判定
 球の方程式– 𝑝は球の中心で𝑟は球の半径。
– 以下を満たすような𝑥は球の表面上の点。

𝑝−𝑥 =𝑟 あるレイとこの球との交差点を計算したいとする– 𝑜はレイの始点で𝑑はレイの方向。
– 𝑃 𝑡 は始点から𝑑の方向に距離𝑡進んだ位置にある点。

𝑃 𝑡 = 𝑜 + 𝑡𝑑

球との交差判定
 球の方程式にレイ上の点の式を代入する– 球とレイが交差する、ということはレイ上の点が球の方程式を満たし、球の表面上に存在するということ。

𝑜

𝑑

𝑃 𝑡1

𝑝 − 𝑃(𝑡) = 𝑟
𝑝

𝑃 𝑡2

𝑝−𝑃 𝑡

𝑝 − 𝑜 − 𝑡𝑑

2

= 𝑟2
2

= 𝑟2

球との交差判定
𝒕についての二次方程式になる– この式の解がレイと球の交差点に対応する。
𝑑 ∙ 𝑑 𝑡 2 − 2𝑑 ∙ 𝑝 − 𝑜 𝑡 + 𝑝 − 𝑜 ∙ 𝑝 − 𝑜 − 𝑟 2 = 0
𝐴 =𝑑∙𝑑𝐵 = −2𝑑 ∙ 𝑝 − 𝑜
𝐶 = 𝑝 − 𝑜 ∙ 𝑝 − 𝑜 − 𝑟2intersect()
二次方程式の判別式を用いて解の存在を判定– さっきの式から判別式を作れる。
判別式

𝐷 𝐵2 − 4𝐴𝐶== 𝑑∙ 𝑝−𝑜
4
4

2

− 𝑑∙𝑑

𝑝 − 𝑜 ∙ 𝑝 − 𝑜 − 𝑟2

ただし、Ray構造体の向きの長さは1にするという約束があったので 𝑑 ∙ 𝑑 = 1 になる。

判別式が負で、解がない場合は交差しないことを意味する

intersect()判別式を使って解を計算する– 二次方程式の解の公式より
𝐵
𝐷
𝐷
−
±
𝑏
±
−𝐵 ± 𝐷
𝐷
2
4
4
𝑡=
=
=
=𝑏±
2𝐴
𝐴
𝐴
4

𝐴 = 𝑑 ∙ 𝑑 = 1という約束
𝐵= −𝑑 ∙ 𝑝 − 𝑜 = −𝑏2intersect() 得られた解について、0.0以下ならレイの向いている側で交差していない

𝑑𝑜 𝑃 𝑡1𝑝
解が負のケース

𝑃 𝑡2

intersect() 得られた解について、kEPS未満なら– レイの始点が、相手の球の表面と非常に近い場所に存在している。
– この場合は球の表面から離れていく方向にレイが発射されていると考え、この交差を無視する。
– レイトレではレイと物体の交差点を始点にして次のレイを追跡する、ということを行うためこのようなことが起きる。

𝑑

𝑃 𝑡1
𝑜
𝑝
解が0.0に近いケース

𝑃 𝑡2

intersect()得られた二つの解について– 両方ともkEPS未満なら交差点は一つもない。

intersect() 得られた二つの解について– 𝒕𝟏 がkEPSより大きいとき• 𝑡1 < 𝑡2 なので𝑡2 もkEPSより大。
• レイと球の交差点は幾何的には二点ある。
• 物理的には始点から発射された
レイは最初にぶつかった点で反射するので、始点に近いほうの交差点、すなわち𝑡1 が交差点になる。

– 𝒕𝟏 がkEPSより小さいとき• 1つ前のif文を通過しているので
𝑡2 はkEPSより大きい。
• よって𝑡2 が交差点。

intersection.h

Hitpoint構造体
 レイと物体の交差点の幾何学的情報を格納する– レイの始点から交差点までの距離 distance– 交差点における物体法線 normal– 交差点のワールド座標系での位置 position
レイ

normaldistance

position

物体表面

Intersection構造体
幾何学的情報（Hitpoint構造体）に加え、物体のIDも格納するintersect()
 交差点までの距離や交差点の位置、交差点における法線をHitpoint構造体に格納する– 球の表面において、交差点における法線は球の中心から交差点へのベクトルを正規化したものになる。
法線𝑝3.3 シーンデータ

scene.hシーンデータレンダリングするシーンのデータ– Sphere構造体の配列– ヘッダに直接記述

シーンデータ球のみでコーネルボックスを表現
ｙｘ

球で囲まれたこの部分がコーネルボックス

intersect_scene() シーンに存在する各物体とレイとの交差判定関数– 各交差点の内、一番レイの始点に近いものを求める。
（一番近い点で反射するから）– 交差判定はすべての物体に対して行う。（全探索）• 計算量はO(N) (Nは物体の数)
レンダリングにおいて最も重い部分の一つはレイトレによる交差判定、すなわちこの関数である。
eduptのシーンデータはかなり単純なので全探索しても問題になりにくいが、交差判定をする対象が何百万個ものポリゴンになったりすると線形の全探索ではもちろんどうにもならない。
こういった場合はkd-tree、BVH、Octreeなどの空間分割木を使って計算量をO(logN)にするのが定石。

4.パストレーシング まえおき– eduptのメインのレンダリング部分の説明の前に光の物理量、レンダリング方程式、そしてレンダリング方程式を解くアルゴリズムの一つパストレーシングについてより一般的な説明をします。

パストレーシング レンダリング方程式に基づいて光の挙動をシミュ
レーションするアルゴリズムの一つ– レンダリング方程式は光の伝達のモデル。
– パストレーシングを使うと写実的なCGを作れる。
• グローバルイルミネーション（光の相互反射）

 光伝達の物理モデルであるレンダリング方程式– 光の物理量について知らなければならない。

4.1 光の物理量

光の物理量物理モデルに従った光の挙動を計算– 光の物理量を把握する必要がある。

放射分析学（Radiometry）– 測光学（Photometry）

放射束(Flux) 単位時間あたり、ある領域を通過するフォトン数（エネルギー）JΦ= Φ[W]s放射照度(Irradiance) 放射束の面積密度（単位面積あたりの放射束）

W𝑑Φ
B(𝑥) 2 =
m
𝑑𝐴

𝑥

平面角
 平面角（平面における角）– 二つの半直線の間の領域。
– 弧の長さと半径の比が角度(radian)。
– 円は2π(rad)の平面角を持つ。

𝑙
𝜃=
𝑟
𝑙

𝜃
𝑟

立体角
 立体角（空間における角）– 球上の領域。
– 球上の領域の面積と半径の二乗比が立体角(steradians)。
– 球は4π(sr)の立体角を持つ。

𝐴
Ω= 2
𝑟
𝑟

Ω

𝐴

放射輝度(Radiance)放射輝度は放射照度の立体角密度（単位立
体角あたりの放射照度）𝑑𝐸𝜔 𝑥
W
𝐿 𝑥, 𝜔 =[ 2]
𝑑𝜔 m ∙ sr
𝜔

𝑑𝜔

𝐸𝜔

𝑥

𝐸𝜔 は𝜔方向に垂直な微小領域における放射照度

ある点においてある方向に通過するフォトンの量が放射輝度と考えることができる。パストレーシングにおいて一つ一つのレイが運ぶ物理量は放射輝度

放射輝度(Radiance)Lambertのコサイン則より𝑑𝐸𝜔 cos 𝜃 = 𝑑𝐸𝑑𝐸 𝑥𝑑2 Φ 𝑥
𝐿 𝑥, 𝜔 ==𝑑𝜔 cos 𝜃 𝑑𝐴𝑑𝜔 cos 𝜃Lambertのコサイン則 領域𝐴をΦ[W]の光が通過しているとすると、領域Φ
𝐴の放射照度は𝐸 = で、領域𝐴𝜔 = 𝐴 cos 𝜃の放射𝐴Φ
Φ
照度は𝐸𝜔 ==𝐴𝜔
𝐴 cos 𝜃 よって、一般に𝐸𝜔 cos 𝜃 = 𝐸が成り立つ。
Φ

𝜔

𝐴𝜔
𝐴

𝜃

4.2 レンダリング方程式

レンダリング方程式
ある空間における光の伝達を記述した方程式– 光の相互反射、グローバルイルミネーションが考慮された式

𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )※真空を仮定

𝐿𝑖
𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

グローバルイルミネーション
 光源から直接届く光の影響だけでなく、他の物体に一回以上反射して間接的に届く光の影響も計算に含めることでグローバルイルミネーションを考慮できる– レンダリング方程式は間接的に届く光も考慮した式になっている
ため、これをきちんと解けばグローバルイルミネーションも考慮
される。

間接光あり

直接光のみ

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )𝒙から𝝎方向への放射輝度
画像の各ピクセルについてカメラからレイを飛ばしシーンとの交差点を𝑥とすると、この𝐿𝑜 を求めることでカメラが受け取る放射輝度が求まる。＝シーンのレンダリングができる。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )𝒙から𝝎方向への自己発光による放射輝度
位置𝑥が光源上にあった場合のみ、この値は０じゃなくなる。Sphere構造体でいうところのemissionがこれ。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )𝒙におけるBRDF
𝜔′ 方向から届き位置𝑥で𝜔方向に反射される光について、反射される程度を示す。放射輝度を放射照度で微分した量として定義される。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )𝝎′方向から𝒙への放射輝度
周囲から𝑥へ降り注ぐ光。この入射光について位置𝑥を中心とした半球上で積分すること
で、位置𝑥に降り注ぐ光の総量を計算できる。それにBRDFをかけることで位置𝑥に降り注ぎ𝜔′ 方向へと反射していく光𝐿𝑜 が求まる。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )コサイン項
面に垂直に入射する光と面に横から入射する光では、同じ放射輝度値でも前者の方が位置𝑥における影響が強い。（極端な話、真横から入射した光は影響が無い）放射照度と放射輝度の関係式に出てくるcos項と本質的には同じ。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

レンダリング方程式
ある空間における光の伝達を記述した方程式
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )𝒙から𝝎方向へのレイとシーンの交差点を𝒓(𝒙, 𝝎)とすると
𝑳𝒊 𝒙, 𝝎′ = 𝑳𝒐 (𝒓 𝒙, 𝝎′ , −𝝎′)という関係が成り立ち
右の積分は再帰的な形になる。

𝐿𝑖

𝜔′

𝐿𝑜
𝜃

𝜔
𝑥

4.3 モンテカルロ積分

モンテカルロ積分
 𝑳𝒐 𝒙, 𝝎 を計算するには複雑な形の積分を解く必要がある
– 再帰的
– 非連続
– 高次元
 モンテカルロ積分– 被積分関数の値を何らかの形で計算できさえすれば解ける。
• 再帰的に評価– 高次元な積分についても誤差の収束速度が次元数の影響を受けない。
• 次元の呪いを受けないモンテカルロ積分
 積分したい式– 𝑓 𝑥 は𝑥が与えられれば計算できるがこの積分を解くこと自体は難しい→モンテカルロ積分
𝐼=𝑓(𝑥)𝑑𝜇(𝑥)𝐷 モンテカルロ積分による推定器
1
𝐼=
𝑁

𝑁

𝑖=1

𝑓 𝑋𝑖
𝑝𝑑𝑓(𝑋𝑖 )サンプル数𝑁𝑋𝑖 は確率変数で確率密度関数pdfに従ってサンプリングされる（もちろん𝑋𝑖 ∈ 𝐷）

モンテカルロ積分
𝑵を増加させていくと推定値𝑰の真値𝑰に対する誤差は減少していく1– 𝑂( )の速度で減少
𝑁

1
𝐼=
𝑁

𝑁

𝑖=1

𝑓 𝑋𝑖
𝑝𝑑𝑓(𝑋𝑖 )Consistentな推定器（アルゴリズム）
𝑁を無限大に大きくしたとき𝐼 が真値𝐼に収束する。
Unbiasedな推定器（アルゴリズム）
𝐼 の期待値が真値𝐼と等しい。

レンダリング方程式に対するモンテカルロ積分積分したい式𝐼=Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )モンテカルロ積分による推定器
1
𝐼=
𝑁

𝑁

𝑖=1

𝑓𝑟 𝑥, 𝜔, 𝜔𝑖′ 𝐿𝑖 𝑥, 𝜔𝑖′ cos 𝜃
𝑝𝑑𝑓(𝜔𝑖′ )サンプル数𝑁𝜔𝑖′ は確率変数で確率密度関数pdfに従ってサンプリングされる。
このケースでは、𝜔𝑖′ は半球上の方向を示すので、半球上で確率密度関数pdfに従って𝑁個の方向をサンプリングすることになる。

レンダリング方程式に対するモンテカルロ積分

合体・展開1
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +
𝑁

𝑁

𝑖=1

𝑓𝑟 𝑥, 𝜔, 𝜔𝑖′ 𝐿𝑖 𝑥, 𝜔𝑖′ cos 𝜃
𝑝𝑑𝑓(𝜔𝑖′ )𝐿𝑖 𝑥, 𝜔′ = 𝐿𝑜 (𝑟 𝑥, 𝜔′ , −𝜔′)より
𝑟 𝑥, 𝜔′ = 𝑥1 とすると再び𝐿𝑜 についてのモンテカルロ積分になり
1
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +𝑁𝑁

𝑓𝑟 𝑥, 𝜔, 𝜔𝑖′ cos 𝜃

𝑖=1

𝑝𝑑𝑓 𝜔𝑖′1(𝐿𝑒 𝑥1 , −𝜔𝑖′ +
𝑁

𝑁

𝑗=1

𝑓𝑟 𝑥1 , −𝜔𝑖′ , 𝜔𝑗′ 𝐿𝑖 𝑥1 , 𝜔𝑗′ cos 𝜃)
𝑝𝑑𝑓(𝜔𝑗′ )同様に
1
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +𝑁
1(𝐿𝑒 𝑥1 , −𝜔𝑖′ +𝑁(𝐿𝑒 𝑥2 , −𝜔𝑖′

1
+
𝑁

𝑁

𝑘=1

𝑁

𝑗=1

𝑁

𝑖=1

𝑓𝑟 𝑥, 𝜔, 𝜔𝑖′ cos 𝜃
𝑝𝑑𝑓 𝜔𝑖′𝑓𝑟 𝑥1 , −𝜔𝑖′ , 𝜔𝑗′ cos 𝜃
𝑝𝑑𝑓 𝜔𝑗′𝑓𝑟 𝑥2 , −𝜔𝑗′ , 𝜔𝑘′ 𝐿𝑖 𝑥2 , 𝜔𝑘′ cos 𝜃))
𝑝𝑑𝑓(𝜔𝑘′ )4.4 パストレーシングパストレーシング レンダリング方程式に対するモンテカルロ積分– 半球上で𝑁個の方向をサンプリングし、その方向にレ
イトレーシング。
– レイとシーンの交差点一つ一つについて、また𝑁個の方向をサンプリングし、その方向にレイトレーシング、を繰り返す。

反射のたびにレイが𝑵倍になる＝
指数的にレイの数が増える！（爆発）パストレーシング
 反射のたびにサンプリングする方向の個数を1個に制限– 指数の底が１になるのでレイの数が指数的に増加することがなくなる。

 パストレーシングによるモンテカルロ積分– 以下をそのままプログラムにしてやればよい。
– 𝐿𝑜 は確率変数になるので、何回も評価して（パストレーシングして）その平均をとることで真値𝐿𝑜 に収束する。

𝑓𝑟 𝑥, 𝜔, 𝜔𝑖′ 𝐿𝑜 (𝑟 𝑥, 𝜔′ , −𝜔′) cos 𝜃𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +𝑝𝑑𝑓(𝜔𝑖′ )パストレーシング 擬似コード 𝑳𝒐 𝒙, 𝝎 ：位置𝒙から𝝎方向への放射輝度1. 位置𝑥を中心とした半球上で1方向サンプリング。（サンプルは確率密度関数pdfに従うとする）
2. 得られたサンプルを𝜔′とする
3. 位置𝑥から𝜔′へレイトレーシング。交差点を𝑥′とする。
4. 𝜃=位置𝑥における法線と𝜔′方向の成す角5.

return 𝐿𝑒 𝑥, 𝜔 +𝑓𝑟 𝑥,𝜔,𝜔′ 𝐿𝑜 𝑥 ′ ,−𝜔′ cos 𝜃
𝑝𝑑𝑓(𝜔′ )

𝑥′

𝜔′ 𝜃
𝑥

パストレーシング 直感的には、光源から出た光がカメラに到達する経路を
カメラ側から逆方向にレイトレーシングによって追跡して求めている、と考えることができる– レイが物体に衝突するたびに半球上で次の方向をサンプリング– 光源に衝突すると、光源からカメラに至る経路の１つが得られた、ということになる– 衝突回数（＝反射回数）は何回相互反射したか、に相当

𝐿𝑒 𝑥, 𝜔レイトレーシングの方向（光源から出た光の進む方向と逆になる）

4.5 ロシアンルーレット

ロシアンルーレット
再帰の終了条件が無い！（＝無限ループ）– 適当な回数再帰したら終了• 統計的Biasが入る。
– 極端な話、再帰回数を一回に限定すると相互反射が一切考慮されない画像がレンダリングされる。真の画像（相互反射が考慮されたもの）と考慮されない画像の差がこの場合のBias。

– ロシアンルーレット• 統計的にUnbiasedなアルゴリズムになる。

ロシアンルーレット
 再帰を続けるか否かをある確率𝒒に基づいて決める– 𝑢を[0, 1]の乱数とし、 0 < 𝑞 ≤ 1とする。
• 𝑢 < 𝑞なら普通に再帰– ただし再帰によって得られた結果を𝑞で割る• さもなければ処理を打ち切る– 期待値を計算すると、ロシアンルーレット導入以前と同じ値になる。
 𝒒の値は任意にとれるが、再帰によって得られる結果に比例させるように決めると分散が小さくなり誤差が小さくなりやすい– eduptでは物体反射率に基づいて決めている。

パストレーシング 擬似コード（ロシアンルーレット） 𝑳𝒐 𝒙, 𝝎 ：位置𝒙から𝝎方向への放射輝度1. 0から1の範囲の乱数を得て𝒒との大小比較
1.
2.

𝒒以上ならreturn 𝑳𝒆 𝒙, 𝝎 として再帰終了
さもなければそのまま処理を続行

2. 位置𝑥を中心とした半球上で1方向サンプリング。（サンプルは確率密度関数pdfに従うとする）
3. 得られたサンプルを𝜔′とする
4. 位置𝑥から𝜔′へレイトレーシング。交差点を𝑥′とする。
5. 𝜃=位置𝑥における法線と𝜔′方向の成す角6.

return 𝐿𝑒 𝑥, 𝜔

𝑓𝑟 𝑥,𝜔,𝜔′ 𝐿𝑜 𝑥 ′ ,−𝜔′ cos 𝜃 𝟏+∙
𝑝𝑑𝑓 𝜔′𝒒ロシアンルーレットの注意 無限ループに陥る確率は極めて低くなるとはいえ、ロシ
アンルーレットをたまたまうまく通過してしまった場合、再帰の深さが深くなりすぎてスタックオーバーフローすることがある– ある一定以上の深さになったらロシアンルーレットの通過確率を極端に下げる。
• 深くなりすぎる確率は激減するが、それでも０ではない。
• Unbiasedなレンダラにこだわる限り、０にすることはできない。
• eduptはこの方法。

– 非再帰パストレーシングにする。
• うまく工夫すると非再帰版に書き直すことが出来る。

– Unbiasedなレンダラを諦める• 再帰の回数に制限を設ければスタックオーバーフローを防ぐことができる（が、Biasが入る）

5.eduptコード解説（レンダリング編）

5.eduptコード解説（レンダリング編）

 まえおきeduptにおけるメインのレンダリング部分の解説。パストレーシングをベースに、光のシーン内での挙動をシミュ
レーション。

eduptファイル概要

main.cpp
edupt::render()を呼び出すだけ。

constant.hintersection.hmaterial.h

各種定数。
交差判定の結果を格納する構造体IntersectionとHitpoint。
Color型と物体材質について記述するReflectionTypeと屈折率の設定。



ppm.h



radiance.hrandom.h



ray.h



render.hscene.h

sphere.hvec.h

レンダリング結果をppm画像として書き出すための関数save_ppm_file()。
ある方向からの放射輝度を得る関数radiance()。パストレーシングのメイン処理。
乱数生成クラスXorShift。
一つ一つの光線、レイを表現する構造体Ray。
レンダリング画像サイズやサンプル数を受け取り、radiance()を使って各ピクセルの具体的な値を決定する関数
render()。
レンダリングするシーンのデータspheres[]とそのシーンに対する交差判定を行う関数intersect_scene()。

基本的な形状としての球を表現するSphere構造体。
ベクトルを表現する構造体Vec。

main.cpp
edupt::render()を呼んでいるだけ– レンダリング解像度– スーパーサンプリング解像度– サブピクセルあたりのサンプリング回数• を指定するrender.h5.1 カメラ設定

カメラの設定カメラ姿勢の設定
camera_up

camera_dir
camera_position

カメラの設定スクリーンの設定
2*screen_widthscreen_dist
camera_position

2*screen_heightカメラの設定スクリーンを張るベクトルの設定
camera_up

screen_ycamera_position

screen_dist
screen_x

camera_dir×camera_up
screen_center画像バッファの作成
レンダリング結果はColor型の配列に格納– 最終的な画像として出力するときは256階調
LDRイメージだが、計算途中ではHDRイメージとしてバッファに保存する。

LDRイメージ
 普通のbmpやpng等の画像はLow Dynamic Range画像– 表現できる輝度の幅が狭い。
– RGBそれぞれ8ビットの階調が一般的。
– 高輝度部分が白飛び。
– 低輝度部分が黒くつぶれる。

白飛び

黒つぶれ

HDRイメージ 物理ベースレンダリングではより表現できる輝度の幅が大きいHigh Dynamic Range画像がしばしば使われる– 物理モデルによる画像は自然界同様非常にダイナミックレンジが広いため。
– RGBそれぞれ8ビットより大きい階調・幅（32ビット等）。
– ディスプレイはLDRによる表示しか対応していないことがほとん
どなので、HDR画像を表示するにはなんらかの方法でLDR画像に変換する必要がある → トーンマッピング

 代表的なフォーマット
– Radiance (.hdr)– OpenEXR (.exr)• http://www.openexr.com/– Floating point TTF

トーンマッピング
eduptではトーンマッピングは深く扱わない– 計算途中ではHDRとして記録しておくも、最終的な保存時には高輝度部分を1.0に切り捨て、輝度値が0.0から1.0の範囲に収まるようにする。

トーンマッピングの例

5.2 画像生成

画像生成

画像生成
 OpenMP– ディレクティブを指定することで簡単に並列化。
– eduptでは画像の行ごとに別々のスレッドで処理させている。

 並列化の処理単位は奥深い問題– なるべく全てのコアを常に100%使用したい。

CGのレンダリングはピクセル間の依存が少なかったり無かったりすることが多い
ため、非常に並列化しやすい問題として知られている。そのため、このように単純なディレクティブを一行指定するだけでもうまく並列化され、高速化が期待で
きる。

画像生成
 画像の各ピクセルについて、走査– 画像の行ごとに並列化される可能性があるため、行ごとに乱数生成機を初期化して使用する。
– 画像の行をまたいで乱数生成機を共有すると競合が発生して結果がおかしくなりうる。

 image_indexがColor配列に対するインデックス

画像生成 スーパーサンプリング– 画像の各ピクセルについて、ピクセルをさらに小領域に分割し、その各領域について輝度値を求める。
– 最終的な各ピクセルの輝度値は小領域全ての平均とする。
サブピクセル

画像上のあるピクセルを
3x3スーパーサンプリングする

サンプリング点（この方向にレイを飛ばして得られた輝度値をこのサブピクセルの色とする）

スーパーサンプリングアンチエイリアスが目的

スーパーサンプリングなし

16x16スーパーサンプリング※ピクセルあたりの総サンプリング回数は同じアンチエイリアス
 ランダムサンプリング

 層化サンプリング（ジッターサンプリング）

 準モンテカルロ法(QMC)
 一様ジッターサンプリング
 etc

画像生成
スクリーン上での各サンプリング点の位置– ピクセルを縦横supersamples個に分割。
– 分割した各小領域＝サブピクセルの中心をサンプリング
点としてその位置を求める。
• その位置からレイを発射、そのピクセルの輝度値を計算すること
になる。

画像生成
 サンプリングする点が決まったらその方向にレイを飛ばし、その方向からの放射輝度をradiance()によって得る– radiance()は確率的な関数なので誤差を含む。
– 同じ方向についてsamples回radiance()を実行し、その平均をその方向からの最終的な放射輝度値とする。

5.3 radiance()radiance() radiance()はレイの方向からの放射輝度値を計算する関数– レイとシーンの交差を𝑥とすると𝐿𝑜 (𝑥, −𝜔)の値を求めればよいことがわかる– 𝐿𝑜 はモンテカルロ積分によって求まるのでradiance()も確率的な関数になるということ
−𝜔 𝑥

𝜔

𝐿𝑜

radiance()
𝐿𝑜 𝑥, 𝜔 = 𝐿𝑒 𝑥, 𝜔 +Ω𝑥

𝑓𝑟 𝑥, 𝜔, 𝜔′ 𝐿𝑖 𝑥, 𝜔′ cos 𝜃 𝑑𝜎(𝜔′ )レンダリング方程式

radiance()のアルゴリズム
 擬似コード（パストレーシングによる𝐿𝑜 の推定を基にしている）
 radiance 𝒙, 𝝎 ：位置𝒙に𝝎方向からくる放射輝度1. 𝑥から𝜔方向にレイトレーシング。交差点を𝑥′とする
2. 0から1の範囲の乱数を得て𝑞との大小比較
1.
2.

3.

4.
5.
6.

𝑞以上ならreturn 𝐿𝑒 𝑥′, 𝜔 として再帰終了
さもなければそのまま処理を続行

位置𝑥′を中心とした半球上で1方向サンプリング。（サンプルは確率密度関数pdfに従うとする）
得られたサンプルを𝜔′とする
𝜃=位置𝑥′における法線と𝜔′方向の成す角
return 𝐿𝑒 𝑥′, 𝜔

𝑓𝑟 𝑥 ′ ,𝜔,𝜔′ 𝑟𝑎𝑑𝑖𝑎𝑛𝑐𝑒(𝑥 ′ ,𝜔′ ) cos 𝜃 1+∙
𝑝𝑑𝑓 𝜔′𝑞radiance.hradiance()
ray方向にレイトレーシングして交差点計算– 交差しなかったら背景と交差したとして
Backgroundcolorを返して終了。

radiance()
orienting_normal– レイの向きと物体法線の内積を計算することで二つのベクトルの成す角のcosθが計算できる。この正負によってレイが物体に外から衝突するのか中から衝突するのかを判定でき、レイに対応する法線orienting_normalを得られる。
hitpoint.normal=orienting_normal
レイ

hitpoint.normal
レイ
orienting_normal

外から衝突するケース

中から衝突するケースradiance()ロシアンルーレット– 無限再帰になるのを防ぐためのロシアンルーレット処理– レイが衝突したオブジェクトの反射率の内最大のものをロシアン
ルーレットの確率にする。
• 反射率が大きければその後の再帰の重要性も高いだろうから処理を打ち切る
確率を下げ、反射率が小さければその後の再帰の重要性も低いだろうから処理を打ち切る確率を上げる。

– スタックオーバーフロー対策でDepthLimit回以上再帰したら打ち
切る確率を上げる。
– いずれにしても、Depth回の再帰は保障する。

物体表面の挙動
 レイが衝突した物体に応じてその後のレイの反射方向を計算する処理が分岐– eduptでは三種類の材質に対応

1. 完全拡散面•

Lambertian

2. 完全鏡面•

いわゆる鏡3. ガラス面•屈折等完全拡散面全方向に等しく反射する面– 入射方向に依らない。
– 拡散反射するときRGBそれぞれ反射する量が違うため物体ごとに異なる色の見え方になる。

完全鏡面
入射角と等しい出射角で反射する面

𝜃 𝜃

𝜃 𝜃

ガラス面
 物体の屈折率の違いによって光が屈折して物体内部に入り込む面– 物体内部に光が入り込むような面では鏡面反射も起きる。
– 屈折方向はSnellの法則で求まる。
– 反射する光の量と屈折する光の量の割合はFresnelの式で求まる。

𝜃 𝜃

屈折率1.0屈折率1.5𝜃′5.4 完全拡散面

radiance()
REFLECTION_TYPE_DIFFUSE– 完全拡散面– 次の反射方向をサンプリングする必要がある。
– orienting_normalの方向を基準として正規直交基底を作り、対応する半球内でサンプリングする。
orienting_normal=wwv

半球内のサンプリング
半球内である方向をサンプリングする– サンプリングに使う確率密度関数は何でも良い。
– 何でも良いが、モンテカルロ積分における被積分関数の形に近い方が分散が小さくなり効率が良くなる。
• インポータンスサンプリング

1
𝐼=
𝑁

𝑁

𝑓 𝑋𝑖
𝑝𝑑𝑓(𝑋𝑖 )𝑖=1
モンテカルロ積分

𝑓𝑟 𝑥 ′ , 𝜔, 𝜔′ 𝑟𝑎𝑑𝑖𝑎𝑛𝑐𝑒(𝑥 ′ , 𝜔′ ) cos 𝜃
被積分関数

完全拡散面におけるインポータンスサンプリング

完全拡散面において
𝜎𝜋– 𝑓𝑟 𝑥 ′ , 𝜔, 𝜔′ 、すなわちBRDFは になる。（𝜎は反射率）– radiance()は未知。
cos 𝜃– 確率密度関数をとして、その確率分布に従って方向𝜋をサンプリングすることにする。

完全拡散面におけるインポータンスサンプリング（ただの計算） 確率密度関数 𝑝𝑑𝑓 𝜃, 𝜙 = 累積分布関数 𝐹 𝜃, 𝜙 =

1
𝜋

cos 𝜃
𝜋

cos 𝜃 𝑑𝜔1 𝜙 𝜃𝐹 𝜃, 𝜙 =cos 𝜃′ sin 𝜃′ 𝑑𝜃′ 𝑑𝜙′
𝜋 0 0
𝜃
1 𝜙
=
𝑑𝜙′ cos 𝜃′ sin 𝜃′ 𝑑𝜃′
𝜋 00𝜃𝜙 cos2 𝜃 ′
=
−
𝜋
2
0

𝜙
=
(1 − cos2 𝜃)
2𝜋
𝐹 𝜙 =

𝜙
, 𝐹 𝜃 = 1 − cos 2 𝜃2𝜋完全拡散面におけるインポータンスサンプリング 確率密度関数 𝑝𝑑𝑓 𝜃, 𝜙 =

cos 𝜃
𝜋

 上のpdfに従って方向をサンプリングする方法– 一つまえのスライドの累積密度関数Fの逆関数を考えると

𝜙𝑖 = 2𝜋𝑢1

𝜃𝑖 = cos −1 𝑢2
ただし𝑢1 , 𝑢2 は0から1の範囲の乱数

 上の式に従って乱数を使ってサンプリング、極座標を使って反射方向dirを計算する。

完全拡散面におけるインポータンスサンプリングcos 𝜃𝜎 確率密度関数がでBRDFが なので最終的にradiance()
𝜋
𝜋

が返すのは

𝐿𝑒 𝑥′, 𝜔 + 𝜎 ∙ 𝑟𝑎𝑑𝑖𝑎𝑛𝑐𝑒(𝑥 ′ , 𝜔′ ) ∙

コードだと

1
𝑞

5.5 完全鏡面・ガラス面

radiance()
REFLECTION_TYPE_SPECULAR– 完全鏡面– 次の反射方向は単純に入射角と等しい出射角の方向
𝐷

𝑁
𝑅

𝜃 𝜃

𝑅 =𝐷−2 𝑁∙𝐷 𝑁

radiance()
REFLECTION_TYPE_REFRACTION– ガラス面– 屈折率は material.hのkIor

 まずSnellの法則によって屈折する方向を求める– 屈折する方向によっては屈折が起こらず、入射した光が全て反射される全反射が起こる。（その場合は完全鏡面と同じ処理になる）

Snellの法則
 入射前の物体の屈折率を𝑛1 、入射後の物体の屈折率を𝑛2 とする。
𝑛1 sin 𝜃1 = 𝑛2 sin 𝜃2 (Snellの法則）

𝑁
𝐷

𝑅

𝜃1 𝜃1

𝑛1
𝑛2

𝜃2

𝑇

𝑛1
2
cos 𝜃2 = 1 −
𝑛2

2

2

(1 − −𝐷 ∙ 𝑁 )右辺が0未満になったとき、𝜃2 が90度を超えたということで、全反射。

Snellの法則
屈折の方向

𝐷

𝑁
𝑅

𝜃1 𝜃1

𝑛1
𝑛2

𝜃2

𝑇

𝑛1
𝑛1
𝑇 = 𝐷−(
𝐷 ∙ 𝑁 + cos2 𝜃2 )𝑁
𝑛2
𝑛2

屈折率
 真空中の光の速度を物質中の光の速度で割った値– 速度が変化することで光は屈折する。
– http://refractiveindex.info/ にたくさん。
水：1.33
水晶：1.54
光学ガラス：1.43-2.14
ダイヤモンド：2.42
サファイア：1.77

 波長によって屈折率は違う– 光の分散Fresnelの式
反射光の運ぶ光の割合はFresnelの式に従う
Schlickの近似が良く使われる– 入射光に対する反射光の運ぶ光の割合𝐹𝑟 𝜃 について入射角度の関数として近似。

𝐷𝐹𝑟 𝜃1 = 𝐹0 + 1 − 𝐹0 1 − cos 𝜃1 5

𝑁

𝜃1 𝜃1

𝑅

𝑛1 > 𝑛2 のときは𝐹𝑟 𝜃1 = 𝐹0 + 1 − 𝐹0 1 − cos 𝜃2 5

𝑛1

𝐹0 は垂直入射における反射の量で𝑛1 − 𝑛2 2
𝐹0 =𝑛1 + 𝑛2 2

𝑛2

𝜃2
𝑇

Fresnelの式
反射光の運ぶ光の割合
𝐹𝑟 𝜃1 = 𝐹0 + 1 − 𝐹0 1 − cos 𝜃1 5
𝑛1 > 𝑛2 のときは𝐹𝑟 𝜃1 = 𝐹0 + 1 − 𝐹0 1 − cos 𝜃2 5

Fresnelの式
屈折光の運ぶ光の割合– 基本的には1 − 𝐹𝑟 𝜃1

レイが運ぶのは放射輝度– 放射輝度は屈折の前後で値が変化する。
– 放射輝度が単位立体角あたりの値なのが理由。
𝑛2 2
屈折直後の放射輝度＝
屈折直前の放射輝度𝑛1ロシアンルーレット
 反射方向と屈折方向の両方から来る放射輝度を計算する
必要がある– 両方に対してradiance()を実行するとレイが指数的に増える→再帰
深度が一定以上ならロシアンルーレットを使って増加を抑制。
– ロシアンルーレットの確率は反射光の運ぶ光の割合に基づいて適当に決める。

通常時
 深度が一定数以下なら普通にレイを追跡– 反射方向、屈折方向のそれぞれから来る放射輝度を計算するためにそれぞれの方向にradiance()を実行。
– Fresnelの式よりすでに計算済みの反射光の運ぶ光の割合と屈折光の運ぶ光の割合をそれぞれに乗算する。

収束について16 sample/pixel
13秒

64 sample/pixel
44秒

256 sample/pixel
174秒

1024 sample/pixel
690秒Intel Core i7 980 (3.33GHz) 10スレッド

4096 sample/pixel
2764秒

参考文献
 Kevin Beason “smallpt http://www.kevinbeason.com/smallpt/”
 David Cline “smallpt presentation
http://www.kevinbeason.com/smallpt/#moreinfo”
 Henrik Wann Jensen (著), 苗村 健 (翻訳) “フォトンマッピング―実写に迫る
コンピュータグラフィックス“ 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”
 Matt Pharr “Image Synthesis (Stanford cs348b) http://candela.stanford.edu/”