レイトレーシング(4): フォトンの生成 #Haskell

光源の定義をやりなおす

前回の記事で、光源について型クラスLightを定義した上で各種光源の型(例えば点光源PointLight)をそのインスタンスとするようにした。その後、シーン情報内で光源のリストを作ろうとした時に問題に気づいた。Haskellでは"同じ型"しかリストにできない。新たに面光源を用意したとして、「光源リスト」に点光源と面光源を混ぜることはできない。オブジェクト指向でのスーパークラスの概念とごちゃまぜにしてしまったのだ。。。

改訂した定義は以下の通り。

Light.hs

data Light = PointLight Color Flux Position3

flux :: Light -> Flux
flux (PointLight _ f _) = f

generatePhoton :: Light -> IO Photon
generatePhoton (PointLight c _ p) = do
  :
(あとは同じ)

引数でPointLightをパターンマッチさせるところがHaskellらしい。今後、平面光源や球光源を追加するときは、dataの定義や各関数のパターンを増やしていけばよい。

フォトンの生成

前回の話に立ち返り、メイン処理の前半を詳細化していこう。メイン処理を次に再掲する。

main = do
  photons <- generatePhotons nphoton lgts
  photoncaches <- tracePhotons objs photons
  a <- forM (photoncaches) $ \i -> do
    putStrLn $ show i
  return ()

このphotons <- generatePhotons nphoton lgtsの部分だ。nphotonとlgtsはそれぞれ、追跡したいフォトン数と光源のリスト。次のように定義しよう。

nphoton = 100000 :: Int

pl1 = PointLight (initColor 1 1 1) 2.0 (initPos 0 3.95 2.5) -- 2W
lgts = [pl1]

フォトン数は生成される画像の品質を左右するパラメータだ。とりあえず10万個としよう。点光源1個で試してみる。この点光源の定義は、赤緑青が同じ比率＝白色光(initColor 1 1 1のところ)であり、発光の強さ(光束、ざっくり言えばエネルギー）が2[W]、光源の位置は x,y,z=0, 3.95, 2.5 [m]ということ。ちなみにSI単位系とし、長さは[m]とする。強さが2[W]だとかなり暗いのではと思われがちだが・・・。光の強さというか明るさについては、白熱電球だと100[W]とか60[W]とかがメジャーだ。ちなみに白熱電球の場合、エネルギー変換効率は2.3〜2.6%ぐらいだそうで、ほとんどが熱などになってしまうということだ。地球に優しくないとして製造中止になるわけだ。一方、最近LED電球ではlumen(lm)を明るさを示すために使っているが白熱電球の100[W]クラスは大体1520[lm]ぐらいらしい。それで、lumenをエネルギーという意味でワットに換算するとき、波長555[nm]の光では、

  1[W] = 683[lm]

だそうなので、波長の違いをとりあえず無視して、1520 ÷ 683 ≒ 2[W] としたのだ。この換算が正しいかどうかよくわからないが。ただし、波長が違うとこの数字が大きく異なるので、本当は無視するとダメかもしれない。

ではあらためてgeneratePhotonsの定義を考えよう。要は、複数の光源からn個のフォトンが放出される、という状況を作り出したい。Lightの定義で、ある光源から1個のフォトンを放出する関数は定義した(generatePhoton :: Light -> IO Photon)。これを使って全光源から放出される数をnphotonにしたいのだ。そうすると次のことを考えていくことになろう。

各光源から放出されるフォトンはそれぞれ何個になるのか?
その計算にはフォトン一個のエネルギーはどれぐらいになるのか?
それがわかれば各光源のエネルギーをフォトン一個のそれで割れば何個かわかるはず!

ということで、コードにしてみた。

Light.hs

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)

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

ここで、powerがフォトン一個のエネルギー(W)でありnsは各光源が放出すべきフォトン数のリストだ。

ここまでくれば後段はその数だけフォトンの生成を繰り返せば良い。具体的には各光源を放出するフォトンの数だけ並べたリストを作り、そのリストの要素毎にgeneratePhotonを呼び出しているだけだ。これで欲しい数のフォトンが得られた。

生成されたフォトンの品質

さて、これで10万個のフォトンを生成したが、果たしてちゃんとランダムに生成できているのだろうか?放出される偏りはないのだろうか?点光源なので、無作為に球状のどの方向にも偏りなく放出されて欲しい。もし完璧に偏りなければ、フォトン数を無限に近づけることで全フォトンの放出方向のベクトル（大きさは1に正規化されている）を足し合わせれば「0ベクトル」になるだろう。既出のメインルーチンの処理を「途中下車」して、生成されたフォトン群を出力させた。さらにそれを下記のプログラムで読み込み、全ベクトルを足し合わせたものを出力してみた。

main :: IO ()
main = do
  np <- getLine
  pw <- getLine
  let nphoton = read np :: Int
      power = read pw :: Double
  dat <- getContents
  pcs <- forM (lines dat) $ \i -> do
    return $ (read i :: Photon)
  let tv = sumVectors pcs
  putStrLn $ show tv

sumVectors :: [Photon] -> Position3
sumVectors [] = o3
sumVectors (pc:pcs) = getDirVec pc + sumVectors pcs

getDirVec :: Photon -> Position3
getDirVec (_, (_, d)) = d

5回実行し、その結果のベクトルの長さを計算したところ下表の通りとなった。もちろん理想は長さゼロだ。

#photon	min	max	avg	err/photon	time
10000	80.5	129.7	101.3	1.01%	0.62 s
100000	47.2	373.4	216.1	0.22%	5.70 s
1000000	531.8	1232.8	792.3	0.079%	57.18 s

これが統計的にみて「十分に無作為」と言えるのかどうか正直微妙だ。しかし今の所、何か改善できるのか、どうしたらできるのかよくわからないのでとりあえず放置しよう。

一方、生成されるフォトンの波長毎の数は無作為か?測ってみたところ下記のようになった。（フォトン数10万個で5回試行）

light	red	green	blue
白色(2W)	33.3%	33.3%	33.3%
白色(2W)+橙色(1W)	44.4%	33.3%	22.2%

橙色の光源は波長の割合を赤2:緑1:青0とした。こちらは綺麗に無作為な生成ができていると言える。まあ、波長については一様乱数で0-1の値を生成してそれを各色の割合で直接選択しているのだからうまくいって当たり前か。

まとめ

とりあえず、フォトンを生成することはできた。次回は各フォトンを追跡してフォトンマップを作ることにしよう。

ソース v.1.0.0.1