レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も） #Haskell

以前に書いたように、フォトンマッピング法は(1)フォトンマップの生成、(2)レイトレーシングの二段階で画像を生成する手法だ。だから、まずはフォトンマップを作らないと始まらない。今回からそのプログラムを作っていこう。

仕様とメインルーチン

まずは仕様から。「シーン」の情報が与えられているとして、次のようなステップで生成したらいいだろう。「シーン」とは描こうとしている三次元世界のことで、情報とは「ここにこういう球がある」とか「ここにスポットライトを配置して」とかいう設定のこととする。

光源のリストから、追跡したいn個のフォトン（光子）を作り出す。
その各フォトンについて、光源からの経路を追跡して物体との衝突記録（フォトンキャッシュと呼ぶことにする）を作成する。
フォトンキャッシュを標準出力へ書き出す。

大枠はこんなに単純だ（本当かどうか知らない、これは実験だから）。さっそくコードにしてみよう。

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

ほとんど上の仕様そのまんま（本当はここに至るまでに、特に最後のforMのところでいろいろ試行錯誤したのだが、出来上がってみると単純だった）。ここで、lgtsは光源のリスト、objsは物体のリストだ。あとは各関数を詳細化していけばよい。Haskellだとこういうトップダウンでの開発がやりやすいように感じる。もちろん、CでもJavaでもそういうアプローチをしているのだが、何をしたいかの次に「どういう手続きにしたらいいか」を考えないといけないのが「素直じゃない」と感じるところかもしれない。あくまで私見だが。Haskellだとこの次に「generatePhotonsは何をするべきか」を記述していけばいいはず。

光源とフォトン

ちょっと横道に逸れて、シーンの情報をどうするか考えておく。上のlgtsやobjsの実際の定義だ。いまのところ"登場人物"は光源と物体だけだ。レイトレーシングの段階になると視点や投影面など他にも必要だが、それは後で考えよう。ではまず光源だ。

以前の記事で、まずは点光源のみ扱うとした。それはそうなのだが、後々面光源なども扱いたいので、型クラスLightを定義して、それを継承する形でPointLight型を作ろう。フォトンマップを作るために、Lightには何が必要か。一つは、その光源の光の強さを提供すること。フォトンが放出される量を知るためだ。もう一つは放出されるフォトンの生成。フォトンの性質（色とか方向とか）は光源の状態に依存するからだ。

前者の関数をflux、後者をgeneratePhotonとしよう。型クラスは次のようになるだろう。

Light.hs

class Light a where
  flux :: a -> Flux
  generatePhoton :: a -> IO Photon

なにやらよく分からない型が登場している。いろいろ試した後だからサクッと書くが、もちろん試行錯誤を繰り返した結果だ。Flux型は発光強度(放射束(radiant flux)、単位は[W])を表すが、実際にはDoubleの別名、Photon型は光源から放出される一つのフォトンを表す。フォトンをどう定義するか大変悩んだが、最終的には次のような仕様にした。

一つのフォトンは特定の波長の光子とする（つまり白色光などではなく単色光）
光源のどこか表面の一点からどこかの方向へ向かう光線とする
反射屈折したフォトンはその地点を起点としてさらに別の方向へ向かう

特に「特定の波長の」が重要だ。一般にレイトレーシング法では光は白色光など「色付き」で処理するはず。様々な波長が混ざった結果を「光」として扱っている。しかし、「複数の波長が混ざった光が物体に衝突すると、果たして反射光や屈折光はどう計算したらいいのか？」という疑問がわく。光が物体に衝突すると、反射、屈折、吸収のいずれかが起こるが、どれぐらい反射するか、どれぐらい吸収するかなどは反射率などで表される。問題はその率が「波長ごとに違う」ことだ。

例えば赤い物体が赤く見えるのは赤以外の光子が吸収されて反射しないから。では、白色光が赤い物体に衝突したら、次は反射したとして追跡を続けるのか吸収として停止するのか？もちろん、各波長の反射率を平均して反射かどうか決めることはできるし、白色光に各波長の反射率を掛け算して反射光としてもよい。しかし、各フォトンの強さが一定のほうがレイトレーシング時の品質に有利と本に書いてあるので掛け算はNGだ。一方で反射率等を平均したら「プリズムで分光された光は表現できるのか？」「赤い球からの相互拡散反射の効果を表現できるのか？」などが疑問だった。

フォトンを特定の波長にするということは、全体的な「色」は複数のフォトンで表現されるということ。白色なら赤、緑、青の少なくとも3つのフォトンがほぼ同数必要だ。白色光ならフォトン1個ですむところ、3個必要なのだ。ならば、追跡するフォトンの数を増やさないと「色がまばらな」画像になってしまうかも知らない。数が増えると大変だが、今回は実験なのだから計算量や時間にはこだわらず、素直な実装で行こう。どうせ、現実世界は特定波長の光子の集まりで照らされているのだから。

あと、返り値がIO Photonなのは、generatePhoton関数の中で乱数を使うから。乱数を使うということは「純粋」じゃなくなるのでIOをつけないとダメらしい。筆者はこれまでHaskellでどのように乱数を生成したらいいのか、IO型を返す関数を実際のプログラムでどのように扱えばいいのかわかってなかったので今回のプログラムに手を出せなかったのだ。ここにきて、randomRIOの存在とIOの使い所がほんの少し分かったので今回手を出してみた。

フォトンの定義

話が長くなった。結局Photon型は、波長＋放射の起点＋放射方向の情報があればよいだろう。波長はとりあえず赤、緑、青のいずれかとしてWavelength型を定義しよう。

Physics.hs

data Wavelength = Red | Green | Blue deriving (Show, Enum)

Enumクラスに属しておけば、色と番号を相互に変換でき、あとあと使えそうだ。

放射の起点と方向だが、これは要するに「直線」のベクトル表現である。ある点$\boldsymbol{p}$を通り方向が$\boldsymbol{d}$である直線$\boldsymbol{r}$は、$t$を任意の実数とすれば

$
\boldsymbol{r} = \boldsymbol{p} + t \cdot \boldsymbol{d}
$

と表される。これはレイトレーシング法で重要な「光線」そのものだ。これをRay型として定義しておこう。Haskellではタプルを使えばよいだろう。

Geometry.hs

type Ray = (Position3, Direction3)

initRay :: Position3 -> Direction3 -> Ray
initRay p d = (p, d)

target :: Double -> Ray -> Position3
target t (p, d) = p + t *> d

getPos :: Ray -> Position3
getPos = fst

getDir :: Ray -> Direction3
getDir = snd

ということでPhoton型はこれらのタプルで定義できる。合わせてフォトンの衝突記録（フォトンキャッシュ）も定義しておこう。衝突記録には、波長、衝突場所、光子の来た方向があればよい。光子の来た方向は、あとあと画像生成時に放射輝度を計算するのに必要そうなので入れておく。結局、フォトンキャッシュもPhotonと同じ構造だ。Rayの意味するところが違うから別の型PhotonCacheとしておこう。

Optics.hs

type Photon = (Wavelength, Ray)

initPhoton :: Wavelength -> Ray -> Photon
initPhoton l r = (l, r)

type PhotonCache = Photon

色

まだ、型の定義が続く・・・。光源の発する光の「色」を指定したい。柔らかいオレンジがかったルームライトもあれば、緑のスポットライトもあるだろう。これをColor型とする。各波長の比率として表せばよさそうだ。

Physics.hs

data Color = Color Double Double Double

initColor :: Double -> Double -> Double -> Color
initColor r g b
  | mag == 0  = Color (1/3) (1/3) (1/3)
  | otherwise = Color (r'/mag) (g'/mag) (b'/mag)
  where
    r' = clipColor r
    g' = clipColor g
    b' = clipColor b
    mag = r' + g' + b'

clipColor :: Double -> Double
clipColor a                  
  | a < 0     = 0
  | a > 1     = 1
  | otherwise = a

decideWavelength :: Color -> Double -> Wavelength
decideWavelength (Color r g b) p
  | p < r              = Red
  | p < r + g          = Green
  | otherwise          = Blue

比率なので、赤緑青の要素を足して1.0になるようにする。黒色は、色の比率はなんでもよくて、光強度がゼロと考える。だから赤緑青が全部ゼロというColorはNGなのだが、万が一指定されたら「白」にする。あと、decideWavelengthは0から1の間の実数を与えたら対応する波長を返す。フォトンの波長をランダムに決めるときに使うものだ。だから、「比率」で表しておくのだ。（エラー処理が面倒なので、0未満なら0、1以上なら1として処理する）

再び光源

長かった。下ごしらえができたので、PointLightの実装に移ろう。点光源なので、持つべき情報は光色、発光強度、位置でいいだろう。

Light.hs

data PointLight = PointLight Color Flux Position3

instance Light PoingLight where
  flux (PointLight _ f _) = f
  generatePhoton (PointLight c _ p) = do
    theta <- randomRIO (0, pi)
    phi   <- randomRIO (0, pi * 2)
    wl    <- randomRIO (0, 1.0)
    let d = initDirFromAngle theta phi
        r = initRay p (fromJust d)
        w = decideWavelength c wl
    return (w, r)

ここで、initDirFromAngleが初出だが説明しておこう。方向ベクトルを生成する関数の極座標版(?)だ。x, y, zを指定する代わりに2つの角度$\theta$と$\phi$を与える。両方ともradianだ。$\theta$はY軸との角度、$\phi$はX軸からZ軸方向への角度とする。点光源は、その位置を中心にあらゆる方向へ均一にフォトンが放出される。ランダムに放出方向を決めるには角度を乱数で指定したほうが楽そうだったのでこのようにした。x, y, zの3つの乱数を使い、条件に合わない場合は棄却する、という方法でランダムな方向ベクトルを得る方法もあるらしい。ただ、よくわからないので無視。

まとめ

今回はフォトンマップ生成プログラムのメインルーチンを考えてみた。たった数行だ。これを今後肉付けしていこう。そのための下ごしらえとして幾つかの型を定義したので、次はこれらを使ってまずはフォトンをn個作る処理を考えようと思う。

(ソース v.1.0.0.1)