CPUの創りかた(2): decorderとmultiplexer
前回は基本的な論理ゲートを作った。今回作ろうとしているCPU (TD4)は以下の本で 説明されているが、論理回路の本格的なところはROMの実装からだ。
ROMの実装においては、まずアドレスを指定するためにdecorderが必要みたいだ。 また後の章では、信号を統合するのにmultiplexerを用いるよう なので、今回は少し複雑な論理回路としてdecorderとmultiplexerを作ろう。
decorder
ガッコウで習ったはずだが、この本を読むまですっかり忘れていた。 decorderは$n$ bitの入力に対し、それを数字と見て数字に該当する出力のみ ONにするようなものだ。今回作るのは負論理なので一つだけL、残りはHになる。 簡単な2 bitの場合の真理値表は次の通り。
A | B | C0 | C1 | C2 | C3 |
---|---|---|---|---|---|
L | L | L | H | H | H |
H | L | H | L | H | H |
L | H | H | H | L | H |
H | H | H | H | H | L |
これを論理回路にすると・・・面倒くさいので本を見て欲しい。 上記の本のp.130に記載してある。
真理値表をそのままにhaskellのコードにしてみよう。
その前に、前回論理回路の入出力の方をBool
の別名としてBin
と定義した。
実際の値はTrue
とFalse
になるが、少し長くてわかりにくいのでsHI
とsLO
を
定義しておく(state HIGH、state LOWのつもり)。
sHI = True :: Bin sLO = False :: Bin -- version 1 lc_decorder2 :: LogicCircuit lc_decorder2 (False:False:_) = [sLO, sHI, sHI, sHI] lc_decorder2 (True:False:_) = [sHI, sLO, sHI, sHI] lc_decorder2 (False:True:_) = [sHI, sHI, sLO, sHI] lc_decorder2 (True:True:_) = [sHI, sHI, sHI, sLO]
できた。とても簡単だ。上の真理値表をそのままコードに落とすだけだ・・・。 しかしROMを作るときは2 bitではなく4 bitが必要らしい。となると16個の 要素を持つリストの定義が16個並ぶわけだ。邪魔くさくて、かつ見にくくなるだろう。 ということで、真理値表を参考に、同等のリストを生成するコードにしてみよう。
まず入力の$n$ bitを数字に変換し、2 bitなら4個、4 bitなら16個の
Bin
要素を持つリストを、変換した数字のところだけL、それ以外はHに
なるよう作ればよいだろう。
入力を数字に直す関数は次の通り。最初の方を下位ビットとして、 上位ビットから1なら足し合わせて二倍する、を繰り返す。
bin2int :: [Bin] -> Int bin2int [] = 0 bin2int (x:xs) = a + 2 * bin2int xs where a = if x == sHI then 1 else 0
これで入力を数字に換えることができた。今度はその数字を$i$とするときに
$i$番目だけLになるようなBin
のリストを作ろう。
decorder' :: Int -> LogicCircuit decorder' b xs = map (\x -> if x == n then sLO else sHI) [0..mx] where mx = 2^b - 1 n = bin2int (take b xs)
第一引数はビット数、そのあとに数字のもとになるビット列を与える。
たとえば、b=2, xs=[sLO, sHI]
だと、2だ(先頭が下位ビット)。
あとは0から3までの数のリストを基に、入力された数と一致するときにL
それ以外はHに変換すればよい。出力は[sHI, sHI, sLO, sHI]
になるはずだ。
これで任意の$n$ bitのdecorderができた!2 bitの場合は以下の様になるだろう。
-- version 2 lc_decorder2 :: LogicCircuit lc_decorder2 = decorder' 2
プログラムらしくなってきた!
・・・いやいや、これは趣旨が違う・・・
今回のネタは、半田ごて等を使った電子工作ができないから、せめて ソフトウェアで論理回路をシミュレーションしようとしたのだ。上記の コードには論理ゲートがこれっぽっちも入っていないではないか! やり直しだ。。。
先に2 bitの場合の回路図を示した。これをそのまま実装しよう! 前回の論理ゲートを使って。
-- version 3 lc_decorder2 :: LogicCircuit lc_decorder2 (a:b:_) = [y0, y1, y2, y3] where [a', b'] = lc_not [a, b] [y0] = lc_nand [a', b'] [y1] = lc_nand [a, b'] [y2] = lc_nand [a', b] [y3] = lc_nand [a, b]
回路図のままなので説明は不要と思う。これだけでも論理回路を作っている気に なるのが楽しいところだ(自分だけかもしれない)。 本当にちゃんと動いているのかテストしよう。先のversion 2と結果を比較するのだ。 同じ名前にできないので、version 2にはダッシュを付けた。
>>> lc_decorder2 [sLO, sLO] == lc_decorder2' [sLO, sLO] True >>> lc_decorder2 [sHI, sLO] == lc_decorder2' [sHI, sLO] True >>> lc_decorder2 [sLO, sHI] == lc_decorder2' [sLO, sHI] True >>> lc_decorder2 [sHI, sHI] == lc_decorder2' [sHI, sHI] True
これをdoctestで実行してみるとちゃんとテストが通る!嬉しい!
さて問題は4 bitだ。回路図がかなり複雑になるが仕方がない (同じく本のp.129に載っている。ただしICの中身として説明されているため、 余分なものもいろいろ含まれている)。これをそのまま実装する。
lc_decorder4 :: LogicCircuit lc_decorder4 (a:b:c:d:_) = [y0, y1, y2 , y3 , y4 , y5 , y6 , y7 ,y8, y9, y10, y11, y12, y13, y14, y15 ] where [a', b', c', d'] = lc_not [a, b, c, d] [a'_b'] = lc_and [a', b'] [a_b' ] = lc_and [a , b'] [a'_b ] = lc_and [a', b ] [a_b ] = lc_and [a , b ] [c'_d'] = lc_and [c', d'] [c_d' ] = lc_and [c , d'] [c'_d ] = lc_and [c', d ] [c_d ] = lc_and [c , d ] [y0] = lc_nand [a'_b', c'_d'] [y1] = lc_nand [a_b' , c'_d'] [y2] = lc_nand [a'_b , c'_d'] [y3] = lc_nand [a_b , c'_d'] [y4] = lc_nand [a'_b', c_d' ] [y5] = lc_nand [a_b' , c_d' ] [y6] = lc_nand [a'_b , c_d' ] [y7] = lc_nand [a_b , c_d' ] [y8] = lc_nand [a'_b', c'_d] [y9] = lc_nand [a_b' , c'_d] [y10] = lc_nand [a'_b , c'_d] [y11] = lc_nand [a_b , c'_d] [y12] = lc_nand [a'_b', c_d] [y13] = lc_nand [a_b' , c_d] [y14] = lc_nand [a'_b , c_d] [y15] = lc_nand [a_b , c_d] lc_decorder4' :: LogicCircuit lc_decorder4' = decorder' 4
実は真理値表をそのまま書いた方が短くて済む気もするのだが、それはそれ、
論理ゲートの組み合わせで作ることに意味があるので。
なお動作確認のため、こちらもdecorder'
を使ったversionを用意した
(lc_decorder4'
)。幾つかテストしてみたが、問題なさそうだ。
>>> lc_decorder4 [sLO, sLO, sLO, sLO] == lc_decorder4' [sLO, sLO, sLO, sLO] True >>> lc_decorder4 [sHI, sLO, sLO, sLO] == lc_decorder4' [sHI, sLO, sLO, sLO] True >>> lc_decorder4 [sHI, sHI, sHI, sHI] == lc_decorder4' [sHI, sHI, sHI, sHI] True
decorderの出来上がりだ!
multiplexer
multiplexerは、複数の入力の中からどれかを選んで出力するものである。 decorderの時もやったように、まずは真理値表と同じ結果が得られるコードを 考えてみよう。
入力値の構造は次のとおりとする。
- チャンネル指定の情報($n$ bit)。2チャンネルなら1 bit、4チャンネルなら2 bit。
- 入力値。2チャンネルなら2 bit。
なので、2チャンネルなら[A,C0,C1]みたいになる。Aの値によってC0かC1が選ばれる。 コードにすると以下のようになるだろう。
multiplexer' :: Int -> LogicCircuit multiplexer' c xs = [xs'!!n] where b = floor (logBase 2 (fromIntegral c)) n = bin2int $ take b xs xs' = drop b xs lc_multiplexer2ch' :: LogicCircuit lc_multiplexer2ch' = multiplexer' 2 lc_multiplexer4ch' :: LogicCircuit lc_multiplexer4ch' = multiplexer' 4
第一引数がチャンネル数、その後は入力値(リスト)だ。チャンネル数から チャンネル指定の情報を取り出すため$log$でビット数を得る。リストの先頭から その分を切り出し、残りのリストの中から$n$番目を選択して返している。
では2チャンネルmultiplexerを論理ゲートで構成してみよう。4チャンネルの 回路図は例によって本のp.176にある。2チャンネルも類推すれば難しくないだろう。 ググってもよい。これまた、回路図をそのままコードにしてみる。
lc_multiplexer2ch :: LogicCircuit lc_multiplexer2ch (a:y0:y1:_) = lc_or (lc_and [a', y0] ++ lc_and [a, y1]) where [a'] = lc_not [a]
論理回路の出力がBin
のリストなので多少ごちゃごちゃしているが、基本的には
回路図のとおり組み合わせていけば良い、というのがミソである。
先のコードと出力を比較テストしてみる。
>>> lc_multiplexer2ch [sLO, sHI, sLO] == lc_multiplexer2ch' [sLO, sHI, sLO] True >>> lc_multiplexer2ch [sHI, sHI, sLO] == lc_multiplexer2ch' [sHI, sHI, sLO] True >>> lc_multiplexer2ch [sLO, sLO, sHI] == lc_multiplexer2ch' [sLO, sLO, sHI] True >>> lc_multiplexer2ch [sHI, sLO, sHI] == lc_multiplexer2ch' [sHI, sLO, sHI] True
うまくいっているようだ。引き続き4チャンネルに進もう。これも回路図の ままに組み合わせていく。
lc_multiplexer4ch :: LogicCircuit lc_multiplexer4ch (a:b:c0:c1:c2:c3:_) = lc_or (y0 ++ y1 ++ y2 ++ y3) where [a', b'] = lc_not [a, b] y0 = lc_and [c0, a', b'] y1 = lc_and [c1, a, b'] y2 = lc_and [c2, a', b] y3 = lc_and [c3, a, b]
multiplexerはdecorderに比べて回路が簡単なので、コードに落としても 簡単だ。ここでは楽するために3入力ANDを使ってしまったが、まあ市販の ICでも3入力はあるみたいなので許してもらおう。
(今回作ったコードはこちら)
まとめ
今回はdecorderとmultiplexerを論理ゲートを組み合わせて作ってみた。 回路図がわかっていれば、それをそのままコードに落とせばいいので ロジックを悩まなくて良い。
decorderができたので、次はROMを実装するかな。
CPUの創りかた(1): 基本論理回路の定義など
さて今回からは、論理回路をやろう。といっても電子工作をする わけではないので、要するにシミュレータを作るということだ。 もちろん大仰なシミュレータではないが。今回唐突に論理回路を 持ち出したのは次の本の「せい」だ。
2003年発行の本で、当時から気になって買いたいと思っていたが そのままになってしまっていた。数ヶ月前にとうとう買ってしまったが、 当方にはその内容がすこぶる面白かった。マシン語は昔やっていたし、 大学で半導体や論理回路(ANDゲートだの加算器だの)も勉強したが、 間を埋めるもの=CPUの構造はわからないままだった。命令はどう フェッチするのか、インストラクションはどこに規定されているのか、 レジスタに値が入ったり出たりするのはどうなっているのか、などなど。 この本には、非常に簡易なCPUだけれどもその辺が明快に記述されている。
となると次は、実際にここに書いてあるCPUを作りたくなるものだ。 が、電子工作の経験はほとんどなく、また道具もないし・・・。 じゃあせめてその論理回路をパソコン上でシミュレートしたいなあ、 というのが今回のネタである。
この本の中では、TD4というCPUを作っている。発売当初はTD4の製作 ネタのHPがいろいろあったようだが今はリンク切れも幾つかある。 ブームは過ぎたかもしれないが、今一度(ソフトウェアで)製作記と 行こう。
論理回路の定義、クロックとか
論理回路といっても単純なANDゲートなどから加算器とか順序回路とか いろいろある。これをどう表現したらいいか? 電源はともかくクロックとかは? などよくわからないことが幾つか出てくる。
Haskellやモナドやライブラリはまだよくわかっていないので、あとあと もっと良い表現・実装方法が見つかればガラッと変えるかもしれないが、 一旦は次のようなものと決めよう。
- 論理回路は複数の入力と複数の出力を持つ関数のようなものとする。
- 端子の状態は0か1で表現されることが多いが、今回はBool型で代用する。 ただしBinという別名をつけておく。
- 入力と出力はBinのリストで表す。
- クロックは厳密に扱わない。論理回路(の関数)を呼び出すことが すなわちクロックが立ち上がったタイミングであるとする。
- 論理回路内の遅延は「無視」する。
- 「状態」を持つ回路については、その状態も出力の一部として出し、 次にその値を入力することで「状態」を表現する。
最後のところはHaskellではStateモナドというものを使えばいいのかも しれないがよくわからないので今の所はこうしておこう。
基本ゲート
早速コードにしてみよう。まず端子の状態を表すBin型を作る。また 論理回路を関数として定義したいので、その引数と返値の型も決めて おこう。こうすることで、その関数は「論理回路」と言える。
type Bin = Bool type LogicCircuit = [Bin] -> [Bin]
試しにAND回路を定義してみる。2入力に限定してもよかったが、今後
多入力もあり得るし、リスト処理関数and
を使えば一緒なので最初から
多入力とした。また入力がない(=空リスト)の時は便宜的にすべての端子が
Lと考えてL=Falseを返すようにした。
-- AND gate (multiple input) lc_and :: LogicCircuit lc_and [] = [False] lc_and xs = [and xs]
2入力で動作確認してみる。
main :: IO () main = do putStrLn ("[L,L]->" ++ (show $ lc_and [False,False])) putStrLn ("[L,H]->" ++ (show $ lc_and [False,True])) putStrLn ("[H,L]->" ++ (show $ lc_and [True,False])) putStrLn ("[H,H]->" ++ (show $ lc_and [True,True]))
実行結果はこれ。
[L,L]->[False] [L,H]->[False] [H,L]->[False] [H,H]->[True]
同様にすればORゲートも簡単だ。ついでにNOTゲートも。
-- OR gate (multiple input) lc_or :: LogicCircuit lc_or [] = [False] lc_or xs = [or xs] -- NOT gate (multiple input) lc_not :: LogicCircuit lc_not [] = [True] lc_not xs = map (not) xs
NOTは多入力というより多数のNOTを一括処理できるイメージだが。。。 それぞれ図にすると下のようになる。
ここまで来るとこれらを使って他のゲートも定義可能だ。電子工作で 頻繁に使われるというNAND、ついでにNORも次の通り。
-- NAND gate (multiple input) lc_nand :: LogicCircuit lc_nand = lc_not . lc_and -- NOR gate (multiple input) lc_nor :: LogicCircuit lc_nor = lc_not . lc_or
NANDはANDの結果を反転させるだけだが、lc_andの返値がリストなので
論理演算子のnot
を使わずlc_not
をかますことにした。ああ簡単だ。
図にすると下の通り。
まとめ
初回なので、まずどのようなシミュレーションをするかを大枠で決め、 次に簡単な論理ゲートを幾つか作ってみた。次はこれらを組み合わせた もう少し大きな回路に取り組んでみたい。
レイトレーシング(9): フィルタ、乱数などで抗ってみる
前回からかなり間が空いてしまった。ほとんど記事は書いていたのだが、 最後の乱数の比較がなかなかできず停滞してしまった。こういう趣味も、 きちんと時間を確保して取り組みたいものだ。。。
さて前回は曲がりなりにもなんとかフォトンマッピング法で画像を生成できた。 ただ、生成画像には課題が多いことも分かった。もちろん作成したのはとても 「限定した仕様」に基づいたものなので仕方のないところもあるが、可能な 限り改善するよう抗ってみる。
と言ってはみたが、結論から述べると今回の取り組みは"そんなに有効でなかった"。 前回示した画像では、球の影がぼやけてしまったり壁の色が均一でないことを 指摘した。「限定した仕様」では影の部分にはまったくフォトンが届かないので、 フォトンの密度から輝度を求めるフォトンマッピング法では「真っ暗」な場所は 表現できないのと、フォトンが少ないとどうしても均一にならないからだ。 と愚痴を言っても仕方がないので、抗った結果を記す。
円錐フィルタ
フォトンマッピング本の7章最後にフィルタについて言及されている。
ここで言うフィルタは、輝度を推定したい点に近いフォトンの重要度を上げ、遠い フォトンは下げることで、その点に近いほど大きい、遠いほど小さい係数を掛けて 合計する。求めたい輝度に関係が少なそうな遠いフォトンの影響を下げようと いうことだ。
円錐フィルタは一次関数的に重み付けを行うもので、計算式は本に記載されている ものをそのまま使った。$w_{pc}$ は次の通り。
$$w_{pc}=1-\frac{d_p}{kr}$$
$k$の値により、フィルタの効き方が変わってくるが、式から考えると $k$は小さくても1.0まで、無限大にするとフィルタの効果が消える。例の 本には$k$が1.1としてあった。フィルタを導入するにあたり、プログラムの 関係する部分を少し変更した。 (ちなみに全ソースはこちら)
estimateRadiance :: Double -> KdTree Double PhotonInfo -> Intersection -> Radiance estimateRadiance pw pmap (p, n, m) | ps == [] = radiance0 | otherwise = (1.0 / (pi * rmax * rmax)) *> (brdf m rad) where ps = filter (isValidPhoton n) $ kNearest pmap nPhoton $ photonDummy p rs = map (\x -> norm ((photonPos x) - p)) ps rmax = maximum rs rad = sumRadiance1 pw rmax rs ps -- Normal (non filter) sumRadiance1 :: Double -> Double -> [Double] -> [PhotonInfo] -> Radiance sumRadiance1 pw rmax rs ps = foldl (+) radiance0 rads where rads = map (photonInfoToRadiance pw) ps -- Cone filter k_cone :: Double k_cone = 1.1 fac_k :: Double fac_k = 1.0 - 2.0 / (3.0 * k_cone) sumRadiance2 :: Double -> Double -> [Double] -> [PhotonInfo] -> Radiance sumRadiance2 pw rmax rs ps = foldl (+) radiance0 rads where wt = map (waitCone (pw / fac_k) rmax) rs rads = zipWith (photonInfoToRadiance) wt ps waitCone :: Double -> Double -> Double -> Double waitCone pw rmax dp = pw * (1.0 - dp / (k_cone * rmax))
estimateRadiance
の一番下、rad=sumRadiance1 pw rmax rs ps
としている。
このsumRadiance1
をフィルタごとに取り替えられるようにするのだ。
フィルタ無しの場合は単にフォトンの出力を足し合わせるだけ(sumRadiance1
)だ。
円錐フィルタでは各フォトンごとに重みを求め('waitCone')、正規化のため
fac_k
で整えている。$k$が1.0, 1.1, 1.5での結果をフィルタ無しと並べ
比較してみる。
奥の壁の下の方、床の色が滲んでいるところがましになっているようだが、球の 影についてはほとんど改善がない。結局影の部分にはフォトンがないため広い 範囲から少しずつフォトンを集めてくる結果、円錐フィルタの傾きが小さくなり、 結果各フォトンの寄与具合はフィルタ無しとたいして変わらなかったのではないか。 結局、注目している交点の近辺に少しでもフォトンがないと効果的でないという ことか?逆に奥の壁は多少だがフォトンがあるため、より効果がわかりやすかった のかもしれない。
試しに、$k$を1.0未満にするとどうなるかだが、0.8と0.67の場合も並べて おいた。0.8で荒れてきて、0.67ではもはや異次元だ・・・。
遠方フォトンの排除
次に試したのがこれ。影の部分では交点の近隣にフォトンがないため、遠方の フォトンまでかき集めてきて放射輝度推定している。要するに影でないところの フォトンを使って影の色を出そうとしているのだから当然無理がある。であれば、 関係ない遠方のフォトンを使わなければいいのでは、という考えだ。 検証では、輝度推定に200個のフォトンを抽出した上で、推定に使うかどうかは 交点からの距離も条件に追加した。
radius2 :: Double radius2 = 0.1 * 0.1 isValidPhoton :: Position3 -> Direction3 -> PhotonInfo -> Bool isValidPhoton p n pi = n <.> (photonDir pi) > 0 && square (p - photonPos pi) < radius2
radius2
が距離の条件である。計算速度を稼ぐため二乗してある。本プログラム
では長さの単位を[m]としており、サンプルのシーンは一辺が4[m]である。この
状況で、距離の条件を0.5、0.2、0.1と変えて試した結果が下図である。
影についてだけ言えば、円錐フィルタより効果はあるように見える。R=0.2
(20[cm])では、古典的レイトレの影に近い感じが出ている。R=0.1
(10[cm])
だと影はかなりくっきりしているが代わりにまだら模様がひどくて如何ともしがたい。
この方法は一定の効果が見込めるものの、強く効かせすぎると推定に使われる フォトンが少なくなってしまい、推定結果がおかしくなりかねない。
メルセンヌツイスター乱数の使用
まだら模様(フォトンが均等に放射されていないことによるノイズ)を減らすためには、 フォトンを多くするか、フォトンの偏りを減らすかだと思う。フォトンを多くすると 計算時間がかかるし、そもそも少ないフォトンでもできるだけ高品質な画像を得たい。 であれば、フォトンの偏りを減らす方向で何ができるだろう?
フォトンの放射方向は乱数を使ってランダムに決めている。コンピュータで真の 乱数を生成するのは困難なので、とても高度な(?)シミュレーションでもない限り、 普通は擬似乱数を使う。その擬似乱数の「乱数としての質」が気になるのだ。 であればより質が良い乱数を使えばどうか?ということでメルセンヌツイスター なる乱数を使ってみることにした。長いので、以後メルセンヌツイスターを MTと書く。
- http://qiita.com/philopon/items/8f647fc8dafe66b7381b
- http://www.math.sci.hiroshima-u.ac.jp/%7Em-mat/MT/mt.html
cabal
でのインストールは下記のようにすれば良い。
$ cabal install mersenne-random
フォトン放射の際に方向ベクトルを生成するが、その関数で使う乱数ライブラリを 取り替えて比較しよう。まずは標準乱数ライブラリ。
import System.Random (中略) generateRandomDir2 :: IO Direction3 generateRandomDir2 = do x <- randomRIO (-1.0, 1.0) y <- randomRIO (-1.0, 1.0) z <- randomRIO (-1.0, 1.0) let v = initPos x y z len = norm v if len > 1.0 || len == 0.0 then generateRandomDir2 else return $ fromJust $ normalize v
次にMT乱数。
import System.Random.Mersenne as MT (中略) generateRandomDir3 :: IO Direction3 generateRandomDir3 = do x' <- MT.randomIO :: IO Double y' <- MT.randomIO :: IO Double z' <- MT.randomIO :: IO Double let x = x' * 2.0 - 1.0 y = y' * 2.0 - 1.0 z = z' * 2.0 - 1.0 v = initPos x y z len = norm v if len > 1.0 || len == 0.0 then generateRandomDir3 else return $ fromJust $ normalize v
実数で生成した場合、0.0-1.0の範囲になるようなので、わざわざ-1.0-1.0に 変換しないといけない。が、まあ大した手間ではない。
標準乱数とMT乱数について性質を比較してみよう。上記の
generateRandomDir2/3
で方向ベクトルを100000個生成し、その方向ベクトルを
原点からの位置ベクトルと見たてて点をプロットしてみた。
(標準乱数)
(MT乱数)
違いがわかるだろうか? 私にはわからない。。。 標準乱数の方はもっとまだら模様だったり特定の箇所に偏っていたりするのかと "期待"していたがそうはならなかった。
標準random
が使っているのはL'Ecuyer(さん?)のアルゴリズムのようだが、
要するに本件で使う程度であれば、そんなに性質の悪いものではないという
ことなのだろう。乱数の良し悪しで言えばこの乱数は周期が小さいなど
MTに比べ劣るところはあるようだが、10万個程度では違いが出てこないのか。
となるとどちらを使っても一緒、下手にライブラリを追加しなくてもいい、 ということになるが処理時間を測ってみると結構な差が出た。 10万個のフォトンを生成するのにかかった処理時間を5回計測した。 いずれもuser時間、単位は秒だ。結果を下表に示す。
乱数ライブラリ | 1 | 2 | 3 | 4 | 5 | 平均 |
---|---|---|---|---|---|---|
標準乱数 | 3.925 | 3.928 | 3.924 | 3.970 | 4.048 | 3.959 |
MT乱数 | 2.355 | 2.318 | 2.471 | 2.529 | 2.543 | 2.443 |
6割ほどMT乱数の方が速い。フォトンマッピング法ではいろいろなところで乱数を 必要とするので、少しでも速い方が良い。 ということで、今後はMT乱数を使うことにしよう。今回唯一の「成果」かな・・・。
まとめ
今回はよりまともに見える画像を生成するよういろいろ試してみたわけだが、 どれもあまりうまくいかなかった。結論としては、
ということかなと。 となれば、俄然(鏡面、拡散)反射に対応してその効果を確かめたい! のだが、レイトレーシングの回も結構続いたので一旦休憩しよう。 ちょっと他のネタをやって、また本プログラムの拡張に取り組みたいと思う。
レイトレーシング(8): 輝度推定による画像生成
前回、画像生成プログラムの大枠を作成した。今回はフォトンマップから各画素の輝度を 計算して画像生成までやってみよう。
各画素の処理
前回の記事でtraceScreen
とした処理の実際の中身を考える。フォトンマッピング法は
光線追跡法をベースにして、(大雑把に言うと)拡散反射光を求めるのにフォトンマップを
使うやりかただ。光源から届く光の量を、光源への光線追跡により求めるかフォトンマップ
から推定するかの違いだ(と思う)。よって、処理の概略は次のとおり。
- 視線と物体との交点を求める。
- 交点における輝度をフォトンマップから推定する。
1の交点を求めるのは、以前にフォトンマップを作成するときに作ったのでそれをそのまま 流用すれば良い。問題は2の輝度の推定だ。
放射輝度推定
1) 交点に近い方からn個のフォトンを抽出する
kdt
ライブラリはkNearest
という関数でn個の近傍点を抽出することができる。
まさにこれを使えば良い。抽出には中心点(交点)を他のフォトン情報と同じ
型で渡さないといけないので交点からダミーのフォトンを作って渡す。p
は交点。
ps = kNearest pmap nPhoton $ photonDummy p photonDummy :: Position3 -> PhotonInfo photonDummy p = PhotonInfo Red p ex3
2) 推定に不要なフォトンを除く(フォトンの入射方向と交点の法線が反対側とか)
このps
から不要なものを除くが、ここでは交点の平面に対し裏側から入射したフォトンを
除外することにしよう。要するに、交点の法線とフォトンの入射方向の内積が負なら省くと
いうことだ。n
は交点の法線ベクトル。
ps' = filter (isValidPhoton n) ps isValidPhoton :: Direction3 -> PhotonInfo -> Bool isValidPhoton n pi = n <.> (photonDir pi) > 0 photonDir :: PhotonInfo -> Direction3 photonDir (PhotonInfo _ _ d) = d
ここで、「不要なフォトンを除外したあとでn個のフォトンを使って推定すべきでは?」 という指摘があるだろう。その通りだと思うが、ちょっと面倒くさいので今回は無視だ。
3,4) 交点と各フォトンの距離を求める+一番遠いフォトンの距離を選ぶ
これは簡単だ。フォトンから位置を取り出し、交点との差を取って大きさを求めればいい。
rs = map (\x -> norm ((photonPos x) - p)) ps' rmax = maximum rs
haskellだと、リスト処理の関数を使えるので簡単だな。
5,6) 各フォトンによる照度を計算して集計する+交点の輝度を計算する
このうち、このステップで使うのはの部分だ。 は、今の所は完全拡散反射しか扱わないことにしたのでとても簡単だ。
これは、完全拡散反射が微小平面で全ての方向へ均等に反射するため半球全体で積分したもの で割って「輝度」にするためだ。 各フォトンから得られた全照度にこれをかければいいのでの外に出せる。 あとはの集計だ。各フォトンに一個分のエネルギーを掛けて 合計してやればよい。
(Qiitaの方でshocker-0x15さんにご指摘いただきました。ありがとうございました。修正しました。)
rad = foldl (+) radiance0 $ map (photonInfoToRadiance pw) ps' photonInfoToRadiance :: Double -> PhotonInfo -> Radiance photonInfoToRadiance pw (PhotonInfo Red _ _) = Radiance pw 0 0 photonInfoToRadiance pw (PhotonInfo Green _ _) = Radiance 0 pw 0 photonInfoToRadiance pw (PhotonInfo Blue _ _) = Radiance 0 0 pw
各フォトンは周波数情報しか持たないので、ここでは赤緑青それぞれに一個のエネルギーを
かけた「輝度」を返しそれを足し合わせている。radiance0
は初期値で「黒」だ。
得られた集計値にを掛け、最後に推定に使った全フォトンが含まれる円の面積で
割って完了だ。円の面積は先の求めた最大半径(rmax
)を使うのだ。なお
を掛けるところは後々の拡張も考え、brdf
という関数に切り出す。
radiande = (1.0 / (pi * rmax * rmax)) *> (brdf m rad) brdf :: Material -> Radiance -> Radiance brdf m rad = (1.0 / pi2) *> ((diffSpec m) <**> rad) diffSpec :: Material -> Color diffSpec (Material d) = d
やっとできた。
サンプルシーンの描画
プログラムができたので、早速実行してみよう。以前作った二種類のフォトンマップで 画像生成した結果は次のとおり。比較のため、使ったフォトンマップ(を可視化したもの)と 古典的光線追跡法で描画したものも並べてみる。
・・・とても微妙な出来栄えだ・・・。今作っているのは「バージョン1」で、点光源だけ、 拡散反射のみ、フォトン追跡は反射を無視(相互拡散反射は計算しない)という条件がついて いるので、ある程度予想していたが、それにしても。。。ちなみに、この条件では古典的 光線追跡法は(たぶん)「完璧な画像」になるはずだ。いわゆる環境光成分も入れていないし。 だから、理想は一番右の画像にすることだ。
例2の生成画像から、ダメダメな点をいくつかピックアップしてみよう。
- 光源から放射される光にバラツキがあるので綺麗な同心円にならない(A)
- 床面の輝度が光が少ない壁面に滲んでいる(B)
- 球の影がくっきり出ていない、ほとんど影になっていない(C)
- 壁が均一の色になっていない(D)
A)について
有限個のフォトンしか使えないのである程度まばらになるのは仕方がないが、もう少し 均一にしたい。そういえば光源から放射されるフォトンは本当に均一なのか疑問だった。 また、層化サンプリングも有効かもしれない。現状は球状にばらまいているだけだが、 球面をいくつかの部分に分けてその部分ごとにランダムに放射すれば均一性が高まるだろう。 次回検証しよう。フォトンの均一性は後述のD)にも影響してくるはずだ。
B)について
交点に近いフォトンを集めて輝度を推定するため致し方ない結果ではある。これを軽減する 手法について例の本に記載があるが。。。今はバージョン1のため特に奥の壁にはフォトンが 到達しないようになっている。機能拡張をしてフォトン追跡で反射も扱うようにすれば、 奥の壁にもフォトンが届き、自ずとこの問題は目立たなくなるだろう。そもそも床面の フォトンを推定に使うことは「問題」ということはなく、それだけのフォトンが近くに 到達しているのだから、計算に入れるのは「妥当」と考えるべきだろう。 この懸案はバージョン2まで保留する。
C)について
一番ショックなのがこれ。影がほとんどない。フォトンマップはそのアルゴリズム上、 フォトンがないと推定できない。このように「まったくフォトンが届かない場所」は 正しく描画できなくて当たり前かもしれない。相互拡散反射に対応したら緩和されるかも しれないが、一方、例の本にある「円錐フィルタ」を適用してその効果を見てみようと おもう。これも次回以降だ。
D)について
これはフォトンマッピング法における「誤差」成分なので、フォトン数を増やしたり、 輝度推定に使うフォトン数を増やせば低減されるものだ。ただフォトン数を増やすにも 限界があるし輝度推定時にフォトンを増やすと推定に使うフォトンの分布範囲が広がって しまうため余計なフォトンを推定に使ってしまい具合が悪い(Cの問題もその一つと言える)。 上の画像は推定に200個のフォトンを使ったが、下記に50,100,300も示そう。
個数を増やすとにじみはかなり消えてくる一方、影はなんだかわからなくなっている。 逆にn=50にすると影は良くなってきているが壁の滲みはかなりひどい。 推定に使うフォトン数の最適値を考えるのは難しいのでとりあえず保留しよう。 まずは他の部分の改良だ。
まとめ
今回、曲がりなりにも画像生成までこぎつけ、古典的光線追跡法と比較することができた。 思ったほど綺麗な画像にならなかったが、それも改善のしがいがあると捉えよう。 次回は光源からのフォトン生成と推定時のフィルタ導入だ。
(今回のソースはこちらから)
レイトレーシング(7): 光線追跡処理の大枠
忙しさにかまけて更新をサボってしまった。。。気を取り直して再開しよう。 前回まででフォトンマップが出来上がったので、今回からはそれを使った光線追跡の プログラムに取り組む。
処理の大枠
フォトンマップを使う以外は古典的な光線追跡のアルゴリズムである。 一番大枠の処理は次の通り。
- 仮想カメラを用意(定義)する(カメラ位置、スクリーン、解像度など)。
- フォトンマップを読み込む。
- スクリーンの各画素について光線追跡を行い輝度を得る。
- 輝度値を画像ファイルに書き出す。
3が本プログラムのメインとなるが、今回は外堀を埋めよう。1、2、4を定義する。 3も画素をループするところまで。
仮想カメラの用意
1について必要なことを書き出してみる。流儀などがあるだろうが、自分が 作るときには以下の要素を用いている。
- 視線
- 視点
- 視線方向: スクリーンの中心方向
- 垂直方向: スクリーンの「上」方向を決めるための情報
- スクリーン
- フォーカス: スクリーンまでの距離としている
- 解像度(x, y)
なお、スクリーンの中心点は視点から視線方向にフォーカス分の距離を進んだ地点とし、 そこから縦横±1.0の領域と固定している。この正方形を解像度の設定値で細かく 分割してそれぞれに光線追跡するわけだ。なお、スクリーンは視線方向に垂直な平面である。
これらの情報は、ちゃんとしたプログラムにするときには設定ファイルなどから読み込んで いろいろ変更できるようにするのだろうが、とにかく「フォトンマッピング」法の効果を 見たいので決め打ち=ソース中に直書きしてしまう。具体的には前回示したフォトン マップの視覚化と同じ設定を使う(あとで画像を比較できるし)。今回使う設定は以下だ。
eyepos = initPos 0 2 0 -- 視点 eyedir = ez3 -- 視線方向 (0, 0, 1)=Z軸 upper = ey3 -- 上方向を表すベクトル (0, 1, 0)=Y軸 focus = 1.0 :: Double -- フォーカス(スクリーンまでの距離) xres = 256 :: Int -- 解像度(横方向) yres = 256 :: Int -- 解像度(縦方向)
ただ、これだけでは光線追跡処理には不十分なので、いくつかの定義を追加する。 最初の画素の位置や次の画素の場所(縦横方向)などを決めてやらないといけないから あらかじめ計算しておく。
stepx = 2.0 / fromIntegral xres :: Double -- 画素の大きさ(横) stepy = 2.0 / fromIntegral yres :: Double -- 画素の大きさ(縦) eex = ex3 -- スクリーンの横方向を表す方向ベクトル(大きさ1) eey = negate ey3 -- スクリーンの縦方向を表す方向ベクトル(大きさ1) origin = focus *> eyedir + ((-1.0 + 0.5 * stepx) *> eex) - (( 1.0 - 0.5 * stepy) *> eey) -- スクリーンの左上の点(最初の画素)
これらを使えば、位置(X, Y)の画素point
は
point = origin + X * stepx * eex + Y * stepy * eey
で表せるわけだ。
フォトンマップの読み込み
さてほとんど忘れかけていたが、前回作ったフォトンマップのフォーマットを思い出してみよう。
次のような形だった。フォトンキャッシュ情報はPhotonCache
型である。
フォトン数(改行) フォトン1個のエネルギー(改行) 1番目のフォトンキャッシュ情報(改行) 2番目のフォトンキャッシュ情報(改行) : n番目のフォトンキャッシュ情報(改行)
よって読み込みはこの逆をしてやればよい。最終的なデータ構造にはk-d treeを使う。 これは例の本(フォトンマッピング)にk-d treeが良いと書いてあったのでそのまま採用 しただけだ。幸い、haskellにはk-d treeのライブラリがいくつかあるので自作する必要は ない。実はここでちょっとつまづいた。ライブラリによって、処理性能(速度)が段違い なのだ。ここでは最終的に採用した kdtを元に進める。最初使ったものは 遅すぎてとても実用にならなかったので(多分、処理速度が数十倍違う)。インストールは cabalから。
cabal install kdt
kdtで使うため、読み込んだデータを変換してPhotonInfo
型のList
にしよう。
PhotonInfo
型の定義、変換用関数は以下のとおり(PhotonCache
型を使えば
いいのだが、フォトンの入射方向を逆転させる必要があるのでこれはあとで取り組む)。
また、k-d treeのデータ構造へ変換するときに必要になるので、フォトンの位置を取り出す
関数も定義しておく。
data PhotonInfo = PhotonInfo Wavelength Position3 Direction3 deriving (Show, Eq) convertToInfo :: PhotonCache -> PhotonInfo convertToInfo (wl, (rp, rd)) = PhotonInfo wl rp (negate rd) infoToPointList :: PhotonInfo -> [Double] infoToPointList (PhotonInfo _ (Vector3 x y z) _) = [x, y, z]
あとは全フォトン情報を読み出してk-d treeに変換すれば良い。なお、フォトンマップの データは標準入力で与えることにする。
各画素の輝度値を求める
具体的な処理は次回以降に回すとして、今回は大枠なので適当な輝度値をでっち上げ、
次の処理でファイルに書き出せるようにしよう。入力は各画素の位置、出力は輝度値とする。
輝度値は、各周波数の値を持った型を定義することにしよう。Radiance
とする。
data Radiance = Radiance Double Double Double deriving (Read, Show)
輝度値も光線追跡の最中に頻繁に演算する。なのでVector
とほとんど同じように
加減などの演算関数を定義することになる。が、まずは型だけ。各画素の位置はxとyの
組みで表そう。そうすると、レイトレーシングのメイン処理は大体次のような形に
なるだろう。全画素の組はあらかじめ作っておき、map
で各組を輝度値に
変換する。
scrmap = [(y, x) | y <- [0..(yres - 1)], x <- [0..(xres - 1)]] : image = map (traceScreen ??) srcmap traceScreen :: <<必要な引数>> -> (Int, Int) -> Radiance traceScreen ?? scrmap = ... :
ひとまず、RGB各周波数の輝度値が 0.05 のRadianceのリストを返すようにしておこう。 あとでわかるが、0.05にしておくと灰色になるはずだ。
輝度値の書き出し
輝度値を画像ファイルへ書き出すため、各周波数成分を256段階の整数に変換したい。 次にような関数を用意すれば良いだろう。
clip :: Double clip = 0.1 gamma = 1.0 / 2.2 rgbmax = 255.0 radianceToRgb :: Double -> Int radianceToRgb d = floor (r * rgbmax) where d' = d / clip r = (if d' > 1.0 then 1.0 else d') ** gamma
なお、clip
は輝度値と最大値"255"(=rgbmax)の対応付けのための閾値。
0.1という設定は、0.1[W]以上のエネルギーを持っていたら最大値とするということ。
この値は適当なので、出来上がる画像の明暗を見て適当に変更するつもりだ。gamma
は
ガンマ補正のための設定値で、一般的(?)な を用いる。
出力する画像データだが、単純なPPM形式を用いよう。フォーマットはググればすぐ わかる。これを標準出力へ吐き出すようにする。以上を踏まえると、書き出し処理の 関数は次のようなものになるだろう。
outputImage :: [Radiance] -> IO () outputImage rs = do putStrLn "P3" putStrLn "## test" putStrLn (show xres ++ " " ++ show yres) putStrLn "255" forM_ rs $ \i -> do putStrLn convertOneCell i
各画素の輝度値のリストを引数にして書きだすだけだ。なおconvertOneCell
は
一画素の結果を文字列に直す関数。
まとめ
上記をまとめ、ソースにしてみた。
(RT.hs) module Main where import System.IO import Control.Monad import Data.KdTree.Static import NumericPrelude import Ray.Algebra import Ray.Geometry import Ray.Physics import Ray.Optics -- PARAMETERS -- -- for camera eyepos = initPos 0 2 0 eyedir = ez3 upper = ey3 focus = 1.0 :: Double xres = 256 :: Int yres = 256 :: Int stepx = 2.0 / fromIntegral xres :: Double stepy = 2.0 / fromIntegral yres :: Double eex = ex3 eey = negate ey3 origin = focus *> eyedir + ((-1.0 + 0.5 * stepx) *> eex) - (( 1.0 - 0.5 * stepy) *> eey) scrmap = [(y, x) | y <- [0..(yres - 1)], x <- [0..(xres - 1)]] -- FUNCTIONS -- main :: IO () main = do (power, photonmap) <- readMap let image = map traceScreen scrmap outputImage image -- readMap :: IO (Double, KdTree Double PhotonInfo) readMap = do np' <- getLine pw' <- getLine let np = read np' :: Int let pw = read pw' :: Double pcs <- forM ([1..np]) $ \i -> do l <- getLine return $ (read l :: PhotonCache) let pmap = build infoToPointList (map convertToInfo pcs) return (pw, pmap) -- traceScreen :: (Int, Int) -> Radiance traceScreen (y, x) = Radiance 0.05 0.05 0.05 -- outputImage :: [Radiance] -> IO () outputImage rs = do putStrLn "P3" putStrLn "## test" putStrLn (show xres ++ " " ++ show yres) putStrLn "255" forM_ rs $ \i -> do putStrLn $ convertOneCell i convertOneCell :: Radiance -> [Char] convertOneCell (Radiance r g b) = (show $ radianceToRgb r) ++ " " ++ (show $ radianceToRgb g) ++ " " ++ (show $ radianceToRgb b)
(Ray/Physics.hs) data Wavelength = Rad | Green | Blue deriving (Show, Read, Enum, Eq)
(Ray/Optics.hs) data Radiance = Radiance Double Double Double deriving (Read, Show) clip = 0.1 :: Double gamma = 1.0 / 2.2 rgbmax = 255.0 radianceToRgb :: Double -> Int radianceToRgb d = floor (r * rgbmax) where d' = d / clip r = (if d' > 1.0 then 1.0 else d') ** gamma type PhotonCache = (Wavelength, Ray) data PhotonInfo = PhotonInfo Wavelength Position3 Direction3 deriving (Show, Eq) infoToPointList :: PhotonInfo -> [Double] infoToPointList (PhotonInfo _ (Vector3 x y z) _) = [x, y, z] convertToInfo :: PhotonCache -> PhotonInfo convertToInfo (wl, (rp, rd)) = PhotonInfo wl rp (negate rd)
コンパイルして実行してみる。
$ ghc -o rt0 RT.hs $ rt0 < photonmap > out1.ppm $ convert out1.ppm out1.png <= ImageMagickを使う
結果は以下のとおり。当たり前だが面白くもなんともない・・・。
ということで、次回は光線追跡の処理を考えていこう。
レイトレーシング(6): やっとフォトン追跡
シーン作成
前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを 作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄に なるので、箱型の空間を用意しよう。これで平面が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
が具体的な
追跡処理なのだが、実行することは次のようになるだろう。
これをコードにしたら次のようになった。
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個ずつ処理を終わらせていけばいいということか。
ような「二重ループ」にすればよいだろう。コードにしてみたのが以下。
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」まで激減した!アルゴリズムは大事だなぁ。(というか最初考えた アルゴリズムがあまりにもイケていなかっただけか・・・)
ようやく完成した。
レイトレーシング(5): 物体の定義
前回までで、光源から放射されるフォトンが生成できたので、次はそれを 追跡してフォトンマップを作ることになるが、そのためには描かれる「物体」を 準備しないといけない。
「物体」定義
次に考える処理は、メインルーチン中の次の部分だ。
photoncaches <- tracePhotons objs photons
objs
はシーン中の物体のリスト、photons
は前回生成したフォトンのリストである。
だから、tracePhotons
の処理を考える前にobjs
がどういうものかを定義する必要がある。
本プログラムでは物体をObject
型で定義しよう。物体に必要な情報は何だろうと考えると
「形」と「材質」だろう、と考えてみた。今後他にも出てくるかもしれないが、今は
その二つにしておく。
data Object = Object Shape Material
材質は木とかガラスとかで、情報としては物体の色や反射率、透明度などなどだ。 考え出すととても複雑な構造になるが、まず作ろうとしている"バージョン1"では 物体は「表面が拡散反射のみ」で「フォトンの追跡は反射を無視」としたのだった。 だからフォトンの追跡は「物体にぶつかったらその位置を記録して終わり」という ことになる。実はこれだとフォトンマップを作るためには材質として何の情報も要らない のだが、あとあとレイトレーシングで画像を作る段になったらさすがに色の情報が 必要なのでその分だけ定義しておこう。「拡散反射率」だ。
data Material = Material Double Double Double
3つのDouble
は赤緑青それぞれの波長での拡散反射率を表す。0から1の間の数値を
設定する。色や輝度をどのような型で表すかまだ決めていないが、その検討次第では
型コンストラクタの引数は型が変わるだろう。
さて、今大事なのは「形」の方だ。レイトレーシングでは、無限平面や球面、二次曲面、 ポリゴンなど様々な「形」が使われる。光源と同じだ。光源の定義では失敗したので、 形の定義では最初から多相性を意識して定義しよう。ただしひとまず無限平面と球面のみ 扱うことにする。(ちょっと先のことも考え、大きさのない"点"も入れておくが)
data Shape = Point Position3 | Plain Direction3 Double | Sphere Position3 Double
ここで、無限平面と球面は次の方程式を満たす三次元空間中の点 (位置ベクトル)の集合である。
ここで、は平面の法線ベクトル、は平面の位置に 関係するパラメータ、は球の中心座標、は球の 半径だ。このあたりの詳しいところはその筋の文献などを参照のこと。たとえば、
- 作者: Fletcher Dunn,Ian Parberry,松田晃一
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/10/04
- メディア: 大型本
- 購入: 21人 クリック: 141回
- この商品を含むブログ (40件) を見る
などに記載がある。上記のPlain
とSphere
の型コンストラクタの引数は、それぞれ
とである。
「形」に必要な関数
フォトンの追跡は、まず物体と衝突する場所(交点)を求めることから始まる。
バージョン1では反射は考えないので交点計算がやることのすべてと言っていい。
Shape型は上記の通り方程式で記述できるので、交点を求めるのは容易だ。光線(Ray
)
との連立方程式を解けばよい。交点を
とすると、
連立させてを求めれば、位置も解るわけだ。
よって、を計算する関数distance
を用意しよう。
distance :: Ray -> Shape -> [Double] -- Point distance r (Point p) = [] -- Plain distance (pos, dir) (Plain n d) | cos == 0 = [] | otherwise = [(d + n <.> pos) / (-cos)] where cos = n <.> dir -- Sphere distance (pos, dir) (Sphere c r) | t1 <= 0.0 = [] | t2 == 0.0 = [t0] | t1 > 0.0 = [t0 - t2, t0 + t2] where o = c - pos t0 = o <.> dir t1 = r * r - (square o - (t0 * t0)) t2 = sqrt t1
交点がない場合、複数の場合があるので、結果はのリストとする。
今後反射や屈折を考えたり、輝度計算をするときには交点での法線ベクトルを求める 必要が出てくる。今はいらないが簡単なので定義しておこう。
getNormal :: Position3 -> Shape -> Maybe Direction3 -- Point getNormal p (Point p') = Nothing -- Plain getNormal p (Plain n d) = Just n -- Sphere getNormal p (Sphere c r) = normalize (p - c)
「点」の場合法線ベクトルがないので、関数getNormal
の結果をMaybe
型に
している。