レイトレーシング(6): やっとフォトン追跡

19628 ワード
PhotonMapping RayTracing Haskell Haskell テキストリンク

シーン作成

前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄になるので、箱型の空間を用意しよう。これで平面が6つ。それだけではつまらないので、1つだけ球面を置こう。その定義を以下に示す。

Scene.hs
m_ball = Material 0.8 0.3 0.3
m_wall = Material 0.8 0.8 0.8
m_ceil = Material 0.4 0.2 0.02
m_flor = Material 0.8 0.6 0.4

wall_bt = initObject (Plain ey3 0) m_flor -- bottom
wall_tp = initObject (Plain (negate ey3) 4) m_ceil
wall_rt = initObject (Plain (negate ex3) 2) m_wall
wall_lt = initObject (Plain ex3 2) m_wall
wall_bk = initObject (Plain ez3 1) m_wall
wall_ft = initObject (Plain (negate ez3) 5) m_wall
ball1   = initObject (Sphere (initPos 1 0.8 3) 0.8) m_ball

objs = [wall_bt, wall_tp, wall_rt, wall_lt, wall_bk, wall_ft, ball1]

横幅と高さが4[m]、奥行き6[m]の部屋に直径1.6[m]の球が置いてあるという想定だ。壁は白、床はオレンジ、球は少し薄めの赤とした。まだ描画していないので色はよくわからないけど。

最終的にobjsが物体のリストである。これを使ってフォトンを追跡しよう。

フォトンの追跡

前回示したtracePhotons関数を定義しよう。といっても、物体のリストとフォトンのリストを引数にするだけだ。

tracePhotons :: [Object] -> [Photon] -> IO [PhotonCache]
tracePhotons objs photons = do
  pcs <- mapM (tracePhoton objs) photons
  return $ concat pcs

うーん、何のことはない。フォトンの数だけtracePhoton(複数形でないことに注意)を呼び出して、結果をフラットなリストにして返しているだけだ。tracePhotonが具体的な追跡処理なのだが、実行することは次のようになるだろう。

  • フォトンをRayとして各物体と衝突する位置=交点までの距離を求める。
  • フォトンの放射された位置から見て前にある交点群について、最も近いものを選ぶ。
  • フォトンキャッシュ情報として返す。

これをコードにしたら次のようになった。

Tracer.hs
tracePhoton :: [Object] -> Photon -> IO [PhotonCache]
tracePhoton os (wl, r) = do
  let iss = filter (\x -> fst x > nearly0) (concat $ map (calcDistance r) os)
      (t, s) = head $ sortBy (comparing fst) iss
  return [(wl, initRay (target t r) (getDir r))]

calcDistance :: Ray -> Object -> [(Double, Object)]
calcDistance r o@(Object s m) = zip ts (replicate (length ts) o)
  where
    ts = distance r s

しつこいようだが"バージョン1"ではフォトンの反射・屈折は扱わないので、上記コードは若干無駄である。例えばtracePhotonの返値はこの場合1個しかないのでIO [PhotonCache]とリストにする必要はない。しかし今後、反射・屈折に対応したら1個のフォトンの追跡で複数のフォトンキャッシュ情報が得られるようになる（ハズ）だから今のうちにその対応をしておこうということだ。また内側で使われているcalcDistanceについても、今必要なのは交点の位置ベクトルだけなのでtarget t rを呼んで距離を位置ベクトルに変換するだけでよいハズだ。しかし、反射・屈折を扱うには交点での材質の情報が必要だ。また反射・屈折方向を求めるためには衝突した物体の「形」も要る。ということで、あえて衝突した位置までの距離と衝突した物体をセットにして返しているわけだ。

フォトンキャッシュ情報の書き出し

さあ、一連の処理の最後だ。各フォトンを追跡して得られたフォトンキャッシュ情報を標準出力に書き出そう。これをファイルに書き込めばフォトンマップのデータファイルとなる。ファイル形式は、

フォトン数(改行)
フォトン1個のエネルギー(改行)
1番目のフォトンキャッシュ情報(改行)
2番目のフォトンキャッシュ情報(改行)
   : 
n番目のフォトンキャッシュ情報(改行)

とする。以前示したメインルーチンより少しシェイプアップして次のようになった。

Main-pm.hs
  putStrLn $ show nphoton
  putStrLn $ show power
  forM_ (photoncaches) $ \i -> do
    putStrLn $ show i

やっとできた!これまでのプログラムをまとめ、コンパイルして実行だ。結果はかなり大きなファイルになる。10万個のフォトンで約130MBぐらい。そうか、反射・屈折なしでこれぐらいということは、それを実装すると簡単に1GBに到達する。まあ、ファイル形式の再考はあとにしよう。

実行結果

ここまでくると、フォトンマップがどんな結果になったのかどうしても見てみたい。もちろん、ちゃんと動いているのか気になる所であるし。フォトンキャッシュ情報にある交点座標を、適当な二次元平面に投影させてフォトンの色も踏まえて描き出してみた結果が以下だ。

おお、なんとなくそれっぽいかも。ちゃんと部屋になって、右下に球も見える。黒いところは球の影でフォトンが届かないところだ。

物は試しに、複数の光源に変えて実行してみる。天井の中央ではなく、左右奥に白と橙色の2つの光源だ。ただ、全体の光強度は同じ、フォトン数も同じだ。

なかなかよい。とりあえず、フォトンマップは完成した、ということにしよう。

メモリ使用量の改善

・・・いや、終わらなかった・・・。

フォトン数を100万個にして処理時間を計測中、たまたま眺めていたtopコマンドの表示に驚いた。処理中の最大メモリ使用量が約1.5GBだったのだ!改めて作成したメインルーチンを見てみよう。

main = do
  (power, photons) <- generatePhotons nphoton lgts
  photoncaches <- tracePhotons objs photons
  putStrLn $ show nphoton
  putStrLn $ show power
  forM_ (photoncaches) $ \i -> do
    putStrLn $ show i

generatePhotons :: Int -> [Light] -> IO (Double, [Photon])
generatePhotons nphoton lights = do
  let power = (sum $ map flux lights) / (fromIntegral nphoton)
      ns    = map (calcN power) lights
  photons <- mapM (mapM generatePhoton) (zipWith replicate ns lights)
  return $ (power, concat photons)

フォトン数に応じていくつも巨大なリストを作っている。フォトンを生成するのにまず、n個の光源のリストを生成し、そのリスト全体に処理をかけている。さらにそのフォトンのリストを処理してフォトンキャッシュのリストを生成、最後にリスト全体を表示している。とにかく大きなリストをいくつも作っている。全フォトンを並列に処理しているのでフォトン数が増えるほどメモリを食うのも仕方なしか。

ここで、フォトンを1個だけ追跡する処理を改めて考えてみよう。

  1. 光源から1個のフォトンを生成する
  2. フォトンを追跡してフォトンキャッシュを求める
  3. フォトンキャッシュを出力する

最後にフォトンキャッシュを出力してしまえば、何も情報を保存しておく必要がないので、並列ではなく直列に1個ずつ処理を終わらせていけばいいということか。

  • まずi番目の光源から放出するフォトン数n(i)を計算し、
  • 各光源iについて一個のフォトン追跡をn(i)回処理する

ような「二重ループ」にすればよいだろう。コードにしてみたのが以下。

main :: IO ()
main = do
  putStrLn $ show nphoton
  let power = (sum $ map flux lgts) / (fromIntegral nphoton)
      ns    = map (calcN power) lgts
  putStrLn $ show power
  zipWithM_ outputPhotonCaches ns lgts

calcN :: Double -> Light -> Int
calcN power light = round (flux light / power)

outputPhotonCaches :: Int -> Light -> IO ()
outputPhotonCaches n lgt = mapM_ (outputPhotonCache lgt) [1..n]

outputPhotonCache :: Light -> Int -> IO ()
outputPhotonCache lgt _ =
  generatePhoton lgt >>= tracePhoton objs >>= mapM_ (putStrLn.show)

mainの中で求めているnsが、各光源から放出されるフォトン数n(i)のリスト。最終行のzipWithM_でそれぞれの光源について処理している(外側のループ)。

outputPhotonCahcesではn(i)回フォトンキャッシュを計算して出力するようにしている(内側のループ)。ちなみに2つ目の引数は使わないので捨てている。このように「同じ処理をn回実行する」というのはどう書くのがスマートなんだろうね。

`outputPhotnCache'(単数形)が実際に一個のフォトンを生成から追跡、最後の書き出しまでやっている関数だ。

これをコンパイルして実行してみると、出力結果は同じだが実行時のメモリ使用量が「数十kB」まで激減した!アルゴリズムは大事だなぁ。（というか最初考えたアルゴリズムがあまりにもイケていなかっただけか・・・）

ようやく完成した。

(ソース v.1.0.0.1)