- レイトレーシング(1): バージョン1の定義、ベクトル演算
- レイトレーシング(2):
Algebraモジュールをいじる - レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も）
- レイトレーシング(4): フォトンの生成
- レイトレーシング(5): 物体の定義
- レイトレーシング(6): やっとフォトン追跡
- レイトレーシング(7): 光線追跡処理の大枠
- レイトレーシング(8): 輝度推定による画像生成
- レイトレーシング(9): フィルタ、乱数などで抗ってみる
シーン作成
前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄になるので、箱型の空間を用意しよう。これで平面が6つ。それだけではつまらないので、1つだけ球面を置こう。その定義を以下に示す。
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として各物体と衝突する位置=交点までの距離を求める。 - フォトンの放射された位置から見て前にある交点群について、最も近いものを選ぶ。
- フォトンキャッシュ情報として返す。
これをコードにしたら次のようになった。
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番目のフォトンキャッシュ情報(改行)
とする。以前示したメインルーチンより少しシェイプアップして次のようになった。
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個ずつ処理を終わらせていけばいいということか。
- まず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」まで激減した!アルゴリズムは大事だなぁ。（というか最初考えたアルゴリズムがあまりにもイケていなかっただけか・・・）
ようやく完成した。