💡

レイトレーシング(12): 集光模様!

2022/09/10に公開

Haskell

Ray Tracing

Photon mapping

tech

(この記事はQiitaからの転載です)

レイトレーシング：目次

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}$ は、

\boldsymbol{v_r} = \boldsymbol{e} - 2 \langle \boldsymbol{e}, \boldsymbol{n} \rangle \boldsymbol{n}

である。ただし、 $\boldsymbol{e}, \boldsymbol{n}$ は正規化(=長さが1)されているものとする。コードは次の通り。

Geometry.hs

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$ とするとそれらの間には以下の関係が成り立つという(理屈はよくわかっていない)。

\eta_1 \sin \theta = \eta_2 \sin \phi

この関係式を使うと、（途中の式変形などはよくわからないが）屈折ベクトル
$\boldsymbol{v_t}$ は次の式で求められる。

\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$ は反射のときと同じ。コードは下記のとおり。

Geometry.hs

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の近似式"というのを代わりに使うそうだ。

F_r \approx F_0 + (1 - F_0) {(1 - \cos \theta)}^5

ここで、 $F_0$ は物質の法線方向からの光の反射率である。物質によりその値が異なる。上記近似式とともに、下記のWebページにいくつかサンプルがある。

ここで、上述の鏡面反射ベクトル計算で返される $\cos \theta$ が必要なのだ。コードはこちら。 $F_0$ (Color型)と $\cos \theta$ が引数である。あとは上記の近似式そのままだ。

Physics.hs

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の各パラメータとロシアンルーレットを使って表現していくのだ。解り易いようにフローチャートで表そう。下図に前回（拡散反射面だけ）と今回（鏡面反射・屈折追加）のフローを示そう。前回のフローは今回のフローの一部になっている(赤点線の部分)のだ。

次にコードで比べてみる。前回のコードはこちら。

Tracer.hs

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は前回取り入れた画質向上策である。

一方、今回拡張したコードはこちら。

Tracer.hs


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. レンダリング方程式(?)

さて、いよいよレンダリングの拡張だ。フォトンマップをみる限り、想定したとおり集光模様が描けそうな気がする！

レンダリングは、注目する点（視線レイと物体の交点）にいろいろな方向から届く光を集計して輝度として返せば良い。ただ、言葉では簡単だがこれがなかなか難しい。私はいまだに良い（正確な/物理的に正しい/表現力の高い/…）集計方法がわからない。とはいえわからないなりに作るしかないので、今回採用した集計方法を説明しよう。

まずは前回までの式。

L_o = L_e + d \frac{1}{\pi} (L_l + L_d)

各文字はそれぞれ以下を表す。

$L_o$ : 交点から視線方向への輝度
$L_e$ : 物体の発光輝度
$L_l$ : 光源からの直接光（フォトンマップから推定することも可）
$L_d$ : 間接光（フォトンマップから推定）
$d$ : 拡散度（Material中のdiffuseness）

これだとそんなに複雑ではない。該当するコードは次の通り。

Tracer.hs

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を適用しているのでちょっとわかりにくいが、基本的に上記の式をそのまま適用している。

鏡面反射と屈折が追加されたことで、輝度計算式は次のように拡張した。

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}$ を掛けているのはなぜだっけ？)

Tracer.hs

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に紫を追加してフォトンの種類を増やして追跡しないとだめだろう。ちょっと面倒なので宿題にしておこう。

ex-11.6-2.png

3. まとめ

今回でやっと集光模様までたどり着いた。最初Haskellでフォトンマッピングによるレイトレーシングプログラムに取り掛かった頃は正直自信がなかったが数年がかりでここまできた。やれやれ。まあ、レイトレーシング・ツールとしてはまだまだやることは大量にあるが、一旦ここまで。気が向いたら拡張しようと思う。

※ 今回分のソースはこちら。

GitHubで編集を提案

1. 今回の改善

-1. Material型

-2. 鏡面反射・屈折ベクトル

(反射)

(屈折)

(Schlickの近似式)

-3. フォトン追跡

-4. レンダリング方程式(?)

-5. おまけ：設定ファイルに対応

2. 作例

3. まとめ

Discussion