- レイトレーシング(1): バージョン1の定義、ベクトル演算
- レイトレーシング(2):
Algebraモジュールをいじる - レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も）
- レイトレーシング(4): フォトンの生成
- レイトレーシング(5): 物体の定義
- レイトレーシング(6): やっとフォトン追跡
- レイトレーシング(7): 光線追跡処理の大枠
- レイトレーシング(8): 輝度推定による画像生成
- レイトレーシング(9): フィルタ、乱数などで抗ってみる
光源の定義をやりなおす
前回の記事で、光源について型クラスLightを定義した上で各種光源の型(例えば点光源PointLight)をそのインスタンスとするようにした。その後、シーン情報内で光源のリストを作ろうとした時に問題に気づいた。Haskellでは"同じ型"しかリストにできない。新たに面光源を用意したとして、「光源リスト」に点光源と面光源を混ぜることはできない。オブジェクト指向でのスーパークラスの概念とごちゃまぜにしてしまったのだ。。。
改訂した定義は以下の通り。
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にしたいのだ。そうすると次のことを考えていくことになろう。
- 各光源から放出されるフォトンはそれぞれ何個になるのか?
- その計算にはフォトン一個のエネルギーはどれぐらいになるのか?
- それがわかれば各光源のエネルギーをフォトン一個のそれで割れば何個かわかるはず!
ということで、コードにしてみた。
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の値を生成してそれを各色の割合で直接選択しているのだからうまくいって当たり前か。
まとめ
とりあえず、フォトンを生成することはできた。次回は各フォトンを追跡してフォトンマップを作ることにしよう。