Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい！

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

graphics

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. 鏡面反射・屈折ベクトル

反射・屈折ベクトルの求め方については、いろいろなところに情報があるので ここでは割愛し、数式とコードのみ示すことにしよう。各記号・文字がわかり 易いように、全体を図にするとこんな感じ。

f:id:eijian:20180707022753p:plain

(反射)

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

f:id:eijian:20180707022738p:plain

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

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はそれぞれ拡散反射、鏡面反射、屈折の場合のフォトン追跡の処理だ。

屈折をサポートしたので、最初の目的だった「集光模様」が得られる はずだ！というわけで、ガラス球を配置してフォトンマップを生成してみよう。

f:id:eijian:20180707022748p:plain

左が不透明な拡散反射面の球の場合、右がガラス球だ。左は拡散反射するので 球の形にフォトンが記録されていて下に球の影が見える。一方右はガラス球 なので球面状にフォトンは記録されず、代わりに球の影の真ん中にフォトンが 集中している部分がある。これが集光模様になるのだ！

-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) ( f L_s + (1 - m) ( (1 - f) L_t ) )

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

  • $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

反射屈折方向にトレースするかどうかの判定など多少ややこしい部分はあるが、 基本的には上記式をそのまま再現している。

さて、それでは拡張した今回のプログラムで数種類の物質の球をレンダリング してみよう。

f:id:eijian:20180707022721p:plain

鏡面反射が加わったことで、かなり表現力が高まった気がする。また、ガラス球が 表現できるようになり念願の集光模様が再現できた！

-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. 作例

上記の他にも少し作例を示そう。

f:id:eijian:20180707022728p:plain

左上は天井からの日光を鏡面三角柱で壁に反射している図、右上は例の フォトンマッピング本でなんども出てくる図。これがやりたかった！ 下の二つは球以外の形状をガラスで作ってみたもの。最後に集光模様の 面白い例として光の波長により屈折率が異なるガラス球の画像。屈折率の違いに よりプリズムのように虹色（？）に分かれた集光模様になる。ただし、本プログラムは RGB三種類のフォトンしか使っていないため、青より短波長の紫が表現できない。 そこまでやるにはRGBに紫を追加してフォトンの種類を増やして追跡しないとだめ だろう。ちょっと面倒なので宿題にしておこう。

f:id:eijian:20180707023003p:plain

3. まとめ

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

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