Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい！

CPUの創りかた(1): 基本論理回路の定義など

logic circuit

さて今回からは、論理回路をやろう。といっても電子工作をする わけではないので、要するにシミュレータを作るということだ。 もちろん大仰なシミュレータではないが。今回唐突に論理回路を 持ち出したのは次の本の「せい」だ。

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を一括処理できるイメージだが。。。 それぞれ図にすると下のようになる。

f:id:eijian:20151213135643p:plain

ここまで来るとこれらを使って他のゲートも定義可能だ。電子工作で 頻繁に使われるという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をかますことにした。ああ簡単だ。 図にすると下の通り。

f:id:eijian:20151213135649p:plain

まとめ

初回なので、まずどのようなシミュレーションをするかを大枠で決め、 次に簡単な論理ゲートを幾つか作ってみた。次はこれらを組み合わせた もう少し大きな回路に取り組んでみたい。

レイトレーシング(9): フィルタ、乱数などで抗ってみる

graphics

前回からかなり間が空いてしまった。ほとんど記事は書いていたのだが、 最後の乱数の比較がなかなかできず停滞してしまった。こういう趣味も、 きちんと時間を確保して取り組みたいものだ。。。

さて前回は曲がりなりにもなんとかフォトンマッピング法で画像を生成できた。 ただ、生成画像には課題が多いことも分かった。もちろん作成したのはとても 「限定した仕様」に基づいたものなので仕方のないところもあるが、可能な 限り改善するよう抗ってみる。

と言ってはみたが、結論から述べると今回の取り組みは"そんなに有効でなかった"。 前回示した画像では、球の影がぼやけてしまったり壁の色が均一でないことを 指摘した。「限定した仕様」では影の部分にはまったくフォトンが届かないので、 フォトンの密度から輝度を求めるフォトンマッピング法では「真っ暗」な場所は 表現できないのと、フォトンが少ないとどうしても均一にならないからだ。 と愚痴を言っても仕方がないので、抗った結果を記す。

円錐フィルタ

フォトンマッピング本の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での結果をフィルタ無しと並べ 比較してみる。

f:id:eijian:20151203173139p:plain

奥の壁の下の方、床の色が滲んでいるところがましになっているようだが、球の 影についてはほとんど改善がない。結局影の部分にはフォトンがないため広い 範囲から少しずつフォトンを集めてくる結果、円錐フィルタの傾きが小さくなり、 結果各フォトンの寄与具合はフィルタ無しとたいして変わらなかったのではないか。 結局、注目している交点の近辺に少しでもフォトンがないと効果的でないという ことか？逆に奥の壁は多少だがフォトンがあるため、より効果がわかりやすかった のかもしれない。

試しに、$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と変えて試した結果が下図である。

f:id:eijian:20151203173151p:plain

影についてだけ言えば、円錐フィルタより効果はあるように見える。R=0.2 (20[cm])では、古典的レイトレの影に近い感じが出ている。R=0.1(10[cm]) だと影はかなりくっきりしているが代わりにまだら模様がひどくて如何ともしがたい。

この方法は一定の効果が見込めるものの、強く効かせすぎると推定に使われる フォトンが少なくなってしまい、推定結果がおかしくなりかねない。

メルセンヌツイスター乱数の使用

まだら模様(フォトンが均等に放射されていないことによるノイズ)を減らすためには、 フォトンを多くするか、フォトンの偏りを減らすかだと思う。フォトンを多くすると 計算時間がかかるし、そもそも少ないフォトンでもできるだけ高品質な画像を得たい。 であれば、フォトンの偏りを減らす方向で何ができるだろう？

フォトンの放射方向は乱数を使ってランダムに決めている。コンピュータで真の 乱数を生成するのは困難なので、とても高度な(?)シミュレーションでもない限り、 普通は擬似乱数を使う。その擬似乱数の「乱数としての質」が気になるのだ。 であればより質が良い乱数を使えばどうか?ということでメルセンヌツイスター なる乱数を使ってみることにした。長いので、以後メルセンヌツイスターを MTと書く。

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個生成し、その方向ベクトルを 原点からの位置ベクトルと見たてて点をプロットしてみた。

(標準乱数) f:id:eijian:20151203173214p:plain

(MT乱数) f:id:eijian:20151203173200p:plain

違いがわかるだろうか? 私にはわからない。。。 標準乱数の方はもっとまだら模様だったり特定の箇所に偏っていたりするのかと "期待"していたがそうはならなかった。

標準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): 輝度推定による画像生成

graphics

前回、画像生成プログラムの大枠を作成した。今回はフォトンマップから各画素の輝度を 計算して画像生成までやってみよう。

各画素の処理

前回の記事でtraceScreenとした処理の実際の中身を考える。フォトンマッピング法は 光線追跡法をベースにして、（大雑把に言うと）拡散反射光を求めるのにフォトンマップを 使うやりかただ。光源から届く光の量を、光源への光線追跡により求めるかフォトンマップ から推定するかの違いだ（と思う）。よって、処理の概略は次のとおり。

  1. 視線と物体との交点を求める。
  2. 交点における輝度をフォトンマップから推定する。
    1. 交点に近い方からn個のフォトンを抽出する
    2. 推定に不要なフォトンを除く(フォトンの入射方向と交点の法線が反対側とか)
    3. 交点と各フォトンの距離を求める
    4. 一番遠いフォトンの距離を選ぶ(最大半径、面積を求めるのに使う）
    5. 各フォトンによる照度を計算して集計する
    6. 交点の輝度を計算する

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) 各フォトンによる照度を計算して集計する+交点の輝度を計算する

書籍「フォトンマッピング」には以下の計算式が書いてある。

{ \displaystyle L_r = \frac{1}{\pi r^{2}} \sum_{p=1}^N f_r \cdot \Delta \Phi_p }

このうち、このステップで使うのは f_r \cdot \Delta \Phi_pの部分だ。 f_rは、今の所は完全拡散反射しか扱わないことにしたのでとても簡単だ。

{ \displaystyle f_r = \frac{(拡散反射率)}{\pi} }

これは、完全拡散反射が微小平面で全ての方向へ均等に反射するため半球全体で積分したもの piで割って「輝度」にするためだ。 各フォトンから得られた全照度にこれをかければいいので \sumの外に出せる。 あとは \Delta \Phi_pの集計だ。各フォトンに一個分のエネルギーを掛けて 合計してやればよい。

(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は初期値で「黒」だ。

得られた集計値に f_rを掛け、最後に推定に使った全フォトンが含まれる円の面積で 割って完了だ。円の面積は先の求めた最大半径(rmax)を使うのだ。なお f_rを掛けるところは後々の拡張も考え、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

やっとできた。

サンプルシーンの描画

プログラムができたので、早速実行してみよう。以前作った二種類のフォトンマップで 画像生成した結果は次のとおり。比較のため、使ったフォトンマップ（を可視化したもの)と 古典的光線追跡法で描画したものも並べてみる。

f:id:eijian:20150822180622p:plain

・・・とても微妙な出来栄えだ・・・。今作っているのは「バージョン1」で、点光源だけ、 拡散反射のみ、フォトン追跡は反射を無視（相互拡散反射は計算しない）という条件がついて いるので、ある程度予想していたが、それにしても。。。ちなみに、この条件では古典的 光線追跡法は（たぶん）「完璧な画像」になるはずだ。いわゆる環境光成分も入れていないし。 だから、理想は一番右の画像にすることだ。

例2の生成画像から、ダメダメな点をいくつかピックアップしてみよう。

f:id:eijian:20150822180630p:plain

  • 光源から放射される光にバラツキがあるので綺麗な同心円にならない(A)
  • 床面の輝度が光が少ない壁面に滲んでいる(B)
  • 球の影がくっきり出ていない、ほとんど影になっていない(C)
  • 壁が均一の色になっていない(D)
A)について

有限個のフォトンしか使えないのである程度まばらになるのは仕方がないが、もう少し 均一にしたい。そういえば光源から放射されるフォトンは本当に均一なのか疑問だった。 また、層化サンプリングも有効かもしれない。現状は球状にばらまいているだけだが、 球面をいくつかの部分に分けてその部分ごとにランダムに放射すれば均一性が高まるだろう。 次回検証しよう。フォトンの均一性は後述のD)にも影響してくるはずだ。

B)について

交点に近いフォトンを集めて輝度を推定するため致し方ない結果ではある。これを軽減する 手法について例の本に記載があるが。。。今はバージョン1のため特に奥の壁にはフォトンが 到達しないようになっている。機能拡張をしてフォトン追跡で反射も扱うようにすれば、 奥の壁にもフォトンが届き、自ずとこの問題は目立たなくなるだろう。そもそも床面の フォトンを推定に使うことは「問題」ということはなく、それだけのフォトンが近くに 到達しているのだから、計算に入れるのは「妥当」と考えるべきだろう。 この懸案はバージョン2まで保留する。

C)について

一番ショックなのがこれ。影がほとんどない。フォトンマップはそのアルゴリズム上、 フォトンがないと推定できない。このように「まったくフォトンが届かない場所」は 正しく描画できなくて当たり前かもしれない。相互拡散反射に対応したら緩和されるかも しれないが、一方、例の本にある「円錐フィルタ」を適用してその効果を見てみようと おもう。これも次回以降だ。

D)について

これはフォトンマッピング法における「誤差」成分なので、フォトン数を増やしたり、 輝度推定に使うフォトン数を増やせば低減されるものだ。ただフォトン数を増やすにも 限界があるし輝度推定時にフォトンを増やすと推定に使うフォトンの分布範囲が広がって しまうため余計なフォトンを推定に使ってしまい具合が悪い（Cの問題もその一つと言える)。 上の画像は推定に200個のフォトンを使ったが、下記に50,100,300も示そう。

f:id:eijian:20150822180643p:plain

個数を増やすとにじみはかなり消えてくる一方、影はなんだかわからなくなっている。 逆にn=50にすると影は良くなってきているが壁の滲みはかなりひどい。 推定に使うフォトン数の最適値を考えるのは難しいのでとりあえず保留しよう。 まずは他の部分の改良だ。

まとめ

今回、曲がりなりにも画像生成までこぎつけ、古典的光線追跡法と比較することができた。 思ったほど綺麗な画像にならなかったが、それも改善のしがいがあると捉えよう。 次回は光源からのフォトン生成と推定時のフィルタ導入だ。

（今回のソースはこちらから）

レイトレーシング(7): 光線追跡処理の大枠

graphics

忙しさにかまけて更新をサボってしまった。。。気を取り直して再開しよう。 前回まででフォトンマップが出来上がったので、今回からはそれを使った光線追跡の プログラムに取り組む。

処理の大枠

フォトンマップを使う以外は古典的な光線追跡のアルゴリズムである。 一番大枠の処理は次の通り。

  1. 仮想カメラを用意（定義）する（カメラ位置、スクリーン、解像度など）。
  2. フォトンマップを読み込む。
  3. スクリーンの各画素について光線追跡を行い輝度を得る。
  4. 輝度値を画像ファイルに書き出す。

3が本プログラムのメインとなるが、今回は外堀を埋めよう。1、2、4を定義する。 3も画素をループするところまで。

仮想カメラの用意

1について必要なことを書き出してみる。流儀などがあるだろうが、自分が 作るときには以下の要素を用いている。

  • 視線
    • 視点
    • 視線方向: スクリーンの中心方向
    • 垂直方向: スクリーンの「上」方向を決めるための情報
  • スクリーン
    • フォーカス: スクリーンまでの距離としている
    • 解像度(x, y)

なお、スクリーンの中心点は視点から視線方向にフォーカス分の距離を進んだ地点とし、 そこから縦横±1.0の領域と固定している。この正方形を解像度の設定値で細かく 分割してそれぞれに光線追跡するわけだ。なお、スクリーンは視線方向に垂直な平面である。

f:id:eijian:20150803220609p:plain

これらの情報は、ちゃんとしたプログラムにするときには設定ファイルなどから読み込んで いろいろ変更できるようにするのだろうが、とにかく「フォトンマッピング」法の効果を 見たいので決め打ち＝ソース中に直書きしてしまう。具体的には前回示したフォトン マップの視覚化と同じ設定を使う(あとで画像を比較できるし)。今回使う設定は以下だ。

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は ガンマ補正のための設定値で、一般的(?)な \gamma = 2.2 を用いる。

出力する画像データだが、単純な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を使う

結果は以下のとおり。当たり前だが面白くもなんともない・・・。

f:id:eijian:20150803220615p:plain

ということで、次回は光線追跡の処理を考えていこう。

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

graphics

シーン作成

前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを 作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄に なるので、箱型の空間を用意しよう。これで平面が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に到達する。 まあ、ファイル形式の再考はあとにしよう。

実行結果

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

f:id:eijian:20150523204302p:plain

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

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

f:id:eijian:20150523204303p:plain

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

メモリ使用量の改善

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

フォトン数を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」まで激減した!アルゴリズムは大事だなぁ。（というか最初考えた アルゴリズムがあまりにもイケていなかっただけか・・・）

ようやく完成した。

レイトレーシング(5): 物体の定義

graphics

前回までで、光源から放射されるフォトンが生成できたので、次はそれを 追跡してフォトンマップを作ることになるが、そのためには描かれる「物体」を 準備しないといけない。

「物体」定義

次に考える処理は、メインルーチン中の次の部分だ。

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

ここで、無限平面と球面は次の方程式を満たす三次元空間中の点\boldsymbol x (位置ベクトル)の集合である。

{ 無限平面: \boldsymbol n \cdot \boldsymbol x = d }

{ 球面: ||\boldsymbol c - \boldsymbol x|| = r }

ここで、 \boldsymbol nは平面の法線ベクトル、 dは平面の位置に 関係するパラメータ、 \boldsymbol cは球の中心座標、 rは球の 半径だ。このあたりの詳しいところはその筋の文献などを参照のこと。たとえば、

実例で学ぶゲーム3D数学

実例で学ぶゲーム3D数学

などに記載がある。上記のPlainとSphereの型コンストラクタの引数は、それぞれ \boldsymbol n, dと \boldsymbol c, rである。

「形」に必要な関数

フォトンの追跡は、まず物体と衝突する場所(交点)を求めることから始まる。 バージョン1では反射は考えないので交点計算がやることのすべてと言っていい。 Shape型は上記の通り方程式で記述できるので、交点を求めるのは容易だ。光線(Ray) との連立方程式を解けばよい。交点を \boldsymbol x = \boldsymbol p + t \cdot \boldsymbol dとすると、 連立させて tを求めれば、位置 \boldsymbol xも解るわけだ。 よって、 tを計算する関数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

交点がない場合、複数の場合があるので、結果は tのリストとする。

今後反射や屈折を考えたり、輝度計算をするときには交点での法線ベクトルを求める 必要が出てくる。今はいらないが簡単なので定義しておこう。

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型に している。

まとめ

今回はフォトン追跡の下準備として、「物体」を定義した。次回はフォトンを追跡 することにしよう。

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

graphics

光源の定義をやりなおす

前回の記事で、光源について型クラス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の値を生成してそれを各色の割合で直接選択しているのだから うまくいって当たり前か。

まとめ

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

※ これまでのソースはここ。 ただし、記事より先に進んでいる場合があるので注意。