0. 以前の記事の訂正
第8回の記事で完全拡散反射の場合のBRDF($f_r$)を下記のように示した。
$
\displaystyle f_r = \frac{1}{2\pi}
$
が、間違いだった。shocker-0x15 さんからご指摘をいただきました。ありがとうございます。正しい式は以下。
$
\displaystyle f_r = \frac{1}{\pi}
$
これにより、得られる輝度が全体的にほぼ倍になるが、輝度を画面の色に変換するところで調整してやればいいので大きな問題はなかった。やれやれ。
1. 今回の改善
今回は例トレーシング回の最終回ということで、やっと目標だった集光模様の実現に取り組むことにする。また、ついでに球以外の形状のサポート、設定ファイルの対応もやってみた。
-1. Material型
まず最初にMaterial型の拡張について説明したい。Material型は物体の物性や模様を表すために定義した型だ。前回までは単純な完全拡散面を想定していたためパラメータは少ししかなかった（コード自体には書いていても使っていないものもあった）。鏡面反射・屈折をサポートするにあたりいくつかのパラメータを追加したし、実際に値を使うようにコードを拡張した。下表に従来と今回の対応を示す。
| parameter | 型 | 説明 | 前回版 | 今回 |
|---|---|---|---|---|
| emittance | Radiance | 自己発光強度 | ○ | ○ |
| reflectance | Color | 拡散反射率(≒物体色?) | ○ | ○ |
| specularRefl | Color | 鏡面反射率(金属反射) | - | ○ |
| transmittance | Color | 透過率 | - | - |
| ior | Color | 屈折率 | - | ○ |
| diffuseness | Double | 拡散度※1 | ○ | ○ |
| metalness | Double | 金属性※2 | - | ○ |
| smoothness | Double | 平滑度 | - | - |
※1: 拡散反射と鏡面反射のあり合いを表す。1.0だと完全拡散反射、0.0だと完全鏡面反射となる。鏡面反射方向からの光をどれだけ反射するか。
※2: 金属性と書いたが金属反射と透過の割合を表す。1.0だと完全に不透明で金属反射のみ、0.0だとガラスのような透明物体となる。
なお、まだ対応できていないパラメータもある。transmittanceは透明物質内での光の通過割合で、色ガラスのような物体を表すため、smoothnessは表面のざらつき具合でぼやけたハイライトなどを表すためのパラメータだ。ぜひともこれらも実装したいが、今後の課題とする。
iorの型がColorになっているのは、色（光の波長）毎に屈折率が違う場合があるため。三角プリズムによる光の分散を再現するために必要だ。（とはいえ、現状は赤、緑、青の3波長しか扱っていないため、いわゆる虹色にならないが）
各パラメータがどのように使われるかは後の節で説明する。
-2. 鏡面反射・屈折ベクトル
反射・屈折ベクトルの求め方については、いろいろなところに情報があるのでここでは割愛し、数式とコードのみ示すことにしよう。各記号・文字がわかり易いように、全体を図にするとこんな感じ。
(反射)
入射ベクトル $\boldsymbol{e}$ 、法線ベクトル $\boldsymbol{n}$ とする時、求める反射ベクトル$\boldsymbol{v_r}$ は、
$ \displaystyle
\boldsymbol{v_r} = \boldsymbol{e} - 2 \langle \boldsymbol{e}, \boldsymbol{n} \rangle \boldsymbol{n}
$
である。ただし、$\boldsymbol{e}, \boldsymbol{n}$ は正規化(=長さが1)されているものとする。コードは次の通り。
specularReflection :: Direction3 -> Direction3 -> (Direction3, Double)
specularReflection n e
| v == Nothing = (n, 0.0)
| c < 0.0 = (fromJust v, -c)
| otherwise = (fromJust v, c)
where
c = e <.> n
v = normalize (e - (2.0 * c) *> n)
まあ、コードは上式をそのまま表現しただけである。戻り値は反射ベクトルの他に内積 $\langle \boldsymbol{e}, \boldsymbol{n} \rangle$ すなわち $\cos \theta$( $\theta$ は$\boldsymbol{e}$ と $\boldsymbol{n}$ のなす角) も含めている。$\cos \theta$ はあとで使われることが多いのだ。
(屈折)
絶対屈折率 $\eta_1$ の物質から絶対屈折率 $\eta_2$ の物質に光が入射した場合、入射角 $\theta$ に対し、屈折ベクトルの角度が法線に対して $\phi$ とするとそれらの間には以下の関係が成り立つという(理屈はよくわかっていない)。
$ \displaystyle
\eta_1 \sin \theta = \eta_2 \sin \phi
$
この関係式を使うと、（途中の式変形などはよくわからないが）屈折ベクトル
$\boldsymbol{v_t}$ は次の式で求められる。
$ \displaystyle
\boldsymbol{v_t} = \frac{\eta_1}{\eta_2} \left( \boldsymbol{e} + \left( \cos \theta - \sqrt{{\left( \frac{\eta_2}{\eta_1} \right)}^2 + {\cos}^2 \theta - 1} \right) \boldsymbol{n} \right)
$
ここで $\cos \theta$ は反射のときと同じ。コードは下記のとおり。
specularRefraction :: Double -> Double -> Double -> Direction3 -> Direction3
-> (Direction3, Double)
specularRefraction ior0 ior1 c0 ed n
| r < 0.0 = (o3, 0.0)
| t == Nothing = (o3, 0.0)
| otherwise = (fromJust t, ior')
where
ior' = ior0 / ior1
r = 1.0 / (ior' * ior') + c0 * c0 - 1.0
a = c0 - sqrt r
n' = if ed <.> n > 0.0 then negate n else n
t = normalize (ior' *> (ed + a *> n'))
引数のior0、ior1、c0は数字が紛らわしいがそれぞれ $\eta_1$, $\eta_2$,$\cos \theta$を意味する。edは入射ベクトル $\boldsymbol{e}$ である。where節を少し説明しよう。屈折ベクトルを求める式中に平方根が含まれている。この中身がマイナスになる場合、光は全反射する（＝屈折ベクトルは存在しない）。
これをチェックするため、一旦平方根の中身をrとし、ガードでマイナスの場合は零ベクトルを返している。さて、屈折ベクトルを求めるには入射光が通ってきたのと同じ側に向いた法線ベクトルが必要だが、本プログラムでは交点での法線は常に物体の「外側」に向くよう定義している。つまり球なら表面から外側に法線ベクトルが向く。しかし屈折では物体内部を光が通って物体表面にて反射と屈折が起きる。この場合、反射屈折ベクトルの計算に必要なのは物体の内側に向いた法線ベクトルである。そのためわざわざ $\cos \theta$ の符号をみて $\boldsymbol{n}$ を反転させている(n')。このことは屈折だけでなく反射も同じだが、ガードによる条件分岐で済ましている。あとは上述の数式通りなので特に問題ないだろう。屈折はおまけとして比屈折率($\frac{\eta_1}{\eta_2}$)も戻り値に入れている。
(Schlickの近似式)
鏡面反射・屈折では、物体表面の反射率が入射角により変化する。正確にはFresnelの公式というので求めるのだが、これがなかなかややこしい。おそらく計算負荷も高い。そこで便利な"Schlickの近似式"というのを代わりに使うそうだ。
$ \displaystyle
F_r \approx F_0 + (1 - F_0) {(1 - \cos \theta)}^5
$
ここで、$F_0$ は物質の法線方向からの光の反射率である。物質によりその値が異なる。上記近似式とともに、下記のWebページにいくつかサンプルがある。
- (https://ja.wikipedia.org/wiki/フレネルの式)
- (http://d.hatena.ne.jp/hanecci/20130525/p3)
- (https://yokotakenji.me/log/math/4501/)
ここで、上述の鏡面反射ベクトル計算で返される $\cos \theta$ が必要なのだ。コードはこちら。$F_0$ (Color型)と $\cos \theta$ が引数である。あとは上記の近似式そのままだ。
reflectionIndex :: Color -> Double -> Color
reflectionIndex (Color r g b) c =
Color (r + (1-r) * c') (g + (1-g) * c') (b + (1-b) * c')
where
c' = (1.0 - c) ** 5.0
-3. フォトン追跡
鏡面反射・屈折が追加されたので、フォトン追跡もそれに合わせていろいろ拡張せねばなるまい。これまでは拡散面だけだったので、その反射率によって拡散反射か吸収（つまり追跡終了）しかなかったが、フォトン追跡の条件分岐のバリエーションが増えた。
- 拡散反射か鏡面反射か吸収か
- 反射か屈折か
この場合分けを前述のMaterialの各パラメータとロシアンルーレットを使って表現していくのだ。解り易いようにフローチャートで表そう。下図に前回（拡散反射面だけ）と今回（鏡面反射・屈折追加）のフローを示そう。前回のフローは今回のフローの一部になっている(赤点線の部分)のだ。
次にコードで比べてみる。前回のコードはこちら。
tracePhoton :: [Object] -> Int -> Photon -> IO [PhotonCache]
tracePhoton os l (wl, r) = do
let is = calcIntersection r os
if is == Nothing
then return []
else do
let (p, n, m) = fromJust is
i <- russianRoulette wl [reflectance m]
pcs <- if i > 0
then reflect p n os l wl
else return []
if (useClassicForDirect == False || l > 0) && diffuseness m > 0.0
then return $ ((wl, initRay p (getDir r)) : pcs)
else return pcs
中程のrussianRouletteにて拡散反射率による分岐をしている。関数reflectは拡散反射してさらにフォトンを追跡するものだ。最後のifは前回取り入れた画質向上策である。
一方、今回拡張したコードはこちら。
tracePhoton :: Bool -> Material -> [Object] -> Int -> Photon
-> IO [PhotonCache]
tracePhoton _ _ _ 10 _ = return []
tracePhoton !uc !m0 !os !l !(wl, r)
| is == Nothing = return []
| otherwise = do
let
is' = fromJust is
(p, _, m) = is'
d = diffuseness m
i <- russianRoulette [d]
ref <- if i > 0
then reflectDiff uc m0 os l wl is'
else reflectSpec uc m0 os l (wl, r) is'
if (uc == False || l > 0) && d > 0.0
then return $ ((wl, initRay p $ getDir r) : ref)
else return ref
where
is = calcIntersection r os
reflectDiff :: Bool -> Material -> [Object] -> Int -> Wavelength
-> Intersection -> IO [PhotonCache]
reflectDiff uc m0 os l wl (p, n, m) = do
i <- russianRoulette [selectWavelength wl $ reflectance m]
if i > 0
then do -- diffuse reflection
dr <- diffuseReflection n
tracePhoton uc m0 os (l+1) $ (wl, initRay p dr)
else return [] -- absorption
reflectSpec :: Bool -> Material -> [Object] -> Int -> Photon -> Intersection
-> IO [PhotonCache]
reflectSpec uc m0 os l (wl, (_, ed)) (p, n, m) = do
let
f0 = selectWavelength wl $ specularRefl m
(rdir, cos0) = specularReflection n ed
f' = f0 + (1.0 - f0) * (1.0 - cos0) ** 5.0
j <- russianRoulette [f']
if j > 0
then tracePhoton uc m0 os (l+1) (wl, initRay p rdir)
else do
if (selectWavelength wl $ ior m) == 0.0
then return [] -- non transparency
else reflectTrans uc m0 os l wl ed (p, n, m) cos0
reflectTrans :: Bool -> Material -> [Object] -> Int -> Wavelength -> Direction3
-> Intersection -> Double -> IO [PhotonCache]
reflectTrans uc m0 os l wl ed (p, n, m) c0 = do
let
ior0 = selectWavelength wl $ ior m0
ior1 = selectWavelength wl $ ior m
(tdir, ior') = specularRefraction ior0 ior1 c0 ed n
m0' = if tdir <.> n < 0.0 then m else m_air
tracePhoton uc m0' os (l+1) (wl, initRay p tdir)
一見して前回と比べてかなり複雑になったのがわかるだろう。russianRouletteがいくつか使われており、それが各分岐になっている。reflectDiff, reflectSpec,reflectTransはそれぞれ拡散反射、鏡面反射、屈折の場合のフォトン追跡の処理だ。
屈折をサポートしたので、最初の目的だった「集光模様」が得られるはずだ！というわけで、ガラス球を配置してフォトンマップを生成してみよう。
左が不透明な拡散反射面の球の場合、右がガラス球だ。左は拡散反射するので球の形にフォトンが記録されていて下に球の影が見える。一方右はガラス球なので球面状にフォトンは記録されず、代わりに球の影の真ん中にフォトンが集中している部分がある。これが集光模様になるのだ！
-4. レンダリング方程式(?)
さて、いよいよレンダリングの拡張だ。フォトンマップをみる限り、想定したとおり集光模様が描けそうな気がする！
レンダリングは、注目する点（視線レイと物体の交点）にいろいろな方向から届く光を集計して輝度として返せば良い。ただ、言葉では簡単だがこれがなかなか難しい。私はいまだに良い（正確な/物理的に正しい/表現力の高い/…）集計方法がわからない。とはいえわからないなりに作るしかないので、今回採用した集計方法を説明しよう。
まずは前回までの式。
$ \displaystyle
L_o = L_e + d \frac{1}{\pi} (L_l + L_d)
$
各文字はそれぞれ以下を表す。
- $L_o$ : 交点から視線方向への輝度
- $L_e$ : 物体の発光輝度
- $L_l$ : 光源からの直接光（フォトンマップから推定することも可）
- $L_d$ : 間接光（フォトンマップから推定）
- $d$ : 拡散度（
Material中のdiffuseness）
これだとそんなに複雑ではない。該当するコードは次の通り。
traceRay :: Int -> Double -> KT.KdTree Double PhotonInfo -> [Object]
-> [Light] -> Ray -> IO Radiance
traceRay 10 _ _ _ _ _ = return radiance0
traceRay l pw pmap objs lgts r
| is == Nothing = return radiance0
| otherwise = return (em + di + ii)
where
is = calcIntersection r objs
(p, n, m) = fromJust is
em = sr_half *> emittance m
di = if useClassicForDirect
then brdf m $ foldl (+) radiance0 $ map (getRadianceFromLight objs p n) lgts
else radiance0
ii = estimateRadiance pw pmap (p, n, m)
emが $L_e$ 、diが $L_l$ 、iiが $L_d$に相当する。直接光と間接光は別々にBRDFを適用しているのでちょっとわかりにくいが、基本的に上記の式をそのまま適用している。
鏡面反射と屈折が追加されたことで、輝度計算式は次のように拡張した。
$ \displaystyle
L_o = \frac{1}{2 \pi} L_e + d \frac{1}{\pi} (L_l + L_d) + (1 - d) \lbrace f L_s + (1 - m) \left((1 - f) L_t \right) \rbrace
$
各文字はそれぞれ以下を表す。
- $L_o$ : 交点から視線方向への輝度
- $L_e$ : 物体の発光輝度
- $L_l$ : 光源からの直接光（フォトンマップから推定することも可）
- $L_d$ : 間接光（フォトンマップから推定）
- $L_s$ : 鏡面反射方向の輝度
- $L_t$ : 屈折方向の輝度
- $d$ : 拡散度（
Material中のdiffuseness） - $f$ : 交点の反射率 (Schlickの近似式で計算)
- $m$ : 金属性（
Material中のmetalness）
なお、$d, f, m$ はいずれも $0 \sim 1$の間の値をとる。コードは次の通り。(そういえば物体の発光輝度に $\frac{1}{2 \pi}$ を掛けているのはなぜだっけ？)
traceRay :: Screen -> Material -> Int -> PhotonMap -> [Object] -> [Light] -> Ray -> IO Radiance
traceRay _ _ 10 _ _ _ _ = return radiance0
traceRay !scr !m0 !l !pmap !objs !lgts !r
| is == Nothing = return radiance0
| otherwise = do
si <- if df == 1.0 || f == black
then return radiance0
else traceRay scr m0 (l+1) pmap objs lgts (initRay p rdir)
ti <- if f' == black || ior1 == 0.0
then return radiance0
else do
let
ior0 = averageIor m0
(tdir, ior') = specularRefraction ior0 ior1 cos0 (getDir r) n
m0' = if tdir <.> n < 0.0 then m else m_air
traceRay scr m0' (l+1) pmap objs lgts (initRay p tdir)
return (sr_half *> emittance m +
df *> brdf m (di + ii) +
(1.0 - df) *> (f <**> si + (1.0 - mt) *> f' <**> ti))
where
is = calcIntersection r objs
(p, n, m) = fromJust is
di = if useClassicForDirect scr
then foldl (+) radiance0 $ map (getRadianceFromLight objs p n) lgts
else radiance0
ii = estimateRadiance scr pmap (p, n, m)
(rdir, cos0) = specularReflection n (getDir r)
df = diffuseness m
mt = metalness m
f = reflectionIndex (specularRefl m) cos0
f' = negateColor f -- this means '1 - f'
ior1 = averageIor m
反射屈折方向にトレースするかどうかの判定など多少ややこしい部分はあるが、基本的には上記式をそのまま再現している。
さて、それでは拡張した今回のプログラムで数種類の物質の球をレンダリングしてみよう。
鏡面反射が加わったことで、かなり表現力が高まった気がする。また、ガラス球が表現できるようになり念願の集光模様が再現できた！
-5. おまけ：設定ファイルに対応
実はこれまではレンダリングする画像の情報（スクリーン情報と物体の情報）はソースコードに直書きしていた。流石にシーンを変更するたびにソースを修正してリコンパイルするのは面倒臭く、またさまざまなシーンを描画するために画像情報を残しておきたいことから、設定ファイルから情報を読み込むように改めた。
まずスクリーン情報（カメラの位置や向き、解像度、アンチエリアス要否など）を表すファイル(Screenファイル)の例。フォーマットは（エセ）YAMLとした。
nphoton : 100000
xresolution : 256
yresolution : 256
antialias : yes # yes or no
samplephoton : 100
useclassic : yes # yes or no
estimateradius: 0.3
ambient : [ 0.001, 0.001, 0.001 ]
maxradiance : 0.01
eyeposition : [ 0.0, 2.0, -4.5 ]
targetposition: [ 0.0, 2.0, 0.0 ]
upperdirection: [ 0.0, 1.0, 0.0 ]
focus : 2.7
photonfilter : none # cone, gauss
こちらは単純な変数名と値の組だけなのでなんのことはない。一方シーン情報（物体の形や配置、色などを定義）の例としてガラス球の場合を示そう。ちょっと長いが勘弁願いたい。
light:
- type : parallelogram
color : [ 1.0, 1.0, 1.0 ]
flux : 5.0
position : [ -0.5, 3.99, 2.5 ]
dir1 : [ 1.0, 0.0, 0.0 ]
dir2 : [ 0.0, 0.0, 1.0 ]
material:
- type : solid
name : mwall
emittance : [ 0.0, 0.0, 0.0 ]
reflectance : [ 0.5, 0.5, 0.5 ]
transmittance: [ 0.0, 0.0, 0.0 ]
specularrefl : [ 0.8, 0.8, 0.8 ]
ior : [ 0.0, 0.0, 0.0 ]
diffuseness : 1.0
metalness : 0.0
smoothness : 0.0
- type : solid
name : mwallr
emittance: [ 0.0, 0.0, 0.0 ]
reflectance: [ 0.4, 0.1, 0.1 ]
transmittance: [ 0.0, 0.0, 0.0 ]
specularrefl: [ 0.0, 0.0, 0.0 ]
ior: [ 0.0, 0.0, 0.0 ]
diffuseness: 1.0
metalness: 0.0
smoothness: 0.0
- type : solid
name: mwallb
emittance: [ 0.0, 0.0, 0.0 ]
reflectance: [ 0.1, 0.1, 0.4 ]
transmittance: [ 0.0, 0.0, 0.0 ]
specularrefl: [ 0.0, 0.0, 0.0 ]
ior: [ 0.0, 0.0, 0.0 ]
diffuseness: 1.0
metalness: 0.0
smoothness: 0.0
- type : solid
name: mparal
emittance: [ 0.7958, 0.7958, 0.7958 ]
reflectance: [ 0.0, 0.0, 0.0 ]
transmittance: [ 0.0, 0.0, 0.0 ]
specularrefl: [ 0.0, 0.0, 0.0 ]
ior: [ 0.0, 0.0, 0.0 ]
diffuseness: 0.0
metalness: 0.0
smoothness: 0.0
- type : solid
name : glass
emittance: [ 0.0, 0.0, 0.0 ]
reflectance: [ 0.0, 0.0, 0.0 ]
transmittance: [ 0.0, 0.0, 0.0 ]
specularrefl: [ 0.08, 0.08, 0.08 ]
ior: [ 1.5, 1.5, 1.5 ]
diffuseness: 0.0
metalness: 0.0
smoothness: 0.0
vertex:
- cl01 : [ -0.5, 3.99, 2.5 ]
- cl02 : [ 0.5, 3.99, 2.5 ]
- cl03 : [ -0.5, 3.99, 3.5 ]
object:
- type : plain
name : flooring
normal : [ 0.0, 1.0, 0.0 ]
position: [ 0.0, 0.0, 0.0 ]
material: mwall
- type : plain
name : ceiling
normal : [ 0.0, -1.0, 0.0 ]
position: [ 0.0, 4.0, 0.0 ]
material: mwall
- type : plain
name : rsidewall
normal : [ -1.0, 0.0, 0.0 ]
position: [ 2.0, 0.0, 0.0 ]
material: mwallb
- type : plain
name : lsidewall
normal : [ 1.0, 0.0, 0.0 ]
position: [ -2.0, 0.0, 0.0 ]
material: mwallr
- type : plain
name : backwall
normal : [ 0.0, 0.0, 1.0 ]
position: [ 0.0, 0.0, -6.0 ]
material: mwall
- type : plain
name : frontwall
normal : [ 0.0, 0.0, -1.0 ]
position: [ 0.0, 0.0, 5.0 ]
material: mwall
- type : sphere
name : ball_glass
center : [ 0.0, 0.8, 3.0 ]
radius : 0.8
material: glass
- type : parallelogram
name : ceiling_light
pos1 : cl01
pos2 : cl02
pos3 : cl03
material: mparal
以前、CPU回で使ったParsecライブラリをここでも使ってみた。CPU回でなんとなくパーサの感触をつかんでいたのであまりハマることなく実装できた。Parsecを使ったパーサ作りは、乗ってくるとどんどんパーサを高度化していけるので楽ちんだ。先人に感謝。
実行方法は次の通り(ガラス球の例、ex-11.6)。
$ cabal build
$ dist/build/pm/pm example/screen0.scr example/ex-11.6.scene | dist/build/rt/rt example/screen0.scr example/ex-11.6.scene | convert - ~/tmp/ex-11.6.png
これで ~/tmp/ex-11.6.pngという画像ファイルが出来上がる。
ただ、同じ情報ファイルをフォトンマップ作成(pm)とレイトレーシング(rt)で
繰り返すのは面倒だし画像フォーマット変換(ImageMagickのconvertを使用)も
毎回書くのは邪魔臭いのでシェルスクリプトにすることにした。
$ util/drv.sh example/screen0.scr example/ex-11.6.scene ~/tmp/ex-11.6.png
2. 作例
上記の他にも少し作例を示そう。
左上は天井からの日光を鏡面三角柱で壁に反射している図、右上は例のフォトンマッピング本でなんども出てくる図。これがやりたかった！下の二つは球以外の形状をガラスで作ってみたもの。最後に集光模様の面白い例として光の波長により屈折率が異なるガラス球の画像。屈折率の違いによりプリズムのように虹色（？）に分かれた集光模様になる。ただし、本プログラムはRGB三種類のフォトンしか使っていないため、青より短波長の紫が表現できない。そこまでやるにはRGBに紫を追加してフォトンの種類を増やして追跡しないとだめだろう。ちょっと面倒なので宿題にしておこう。
3. まとめ
今回でやっと集光模様までたどり着いた。最初Haskellでフォトンマッピングによるレイトレーシングプログラムに取り掛かった頃は正直自信がなかったが数年がかりでここまできた。やれやれ。まあ、レイトレーシング・ツールとしてはまだまだやることは大量にあるが、一旦ここまで。気が向いたら拡張しようと思う。
※ 今回分のソースはこちら。