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

(転載その9)

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

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

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

円錐フィルタ

フォトンマッピング本(フォトンマッピング、Henrik Wann Jensen著、オーム社)の7章最後にフィルタについて言及されている。ここで言うフィルタは、輝度を推定したい点に近いフォトンの重要度を上げ、遠いフォトンは下げることで、その点に近いほど大きい、遠いほど小さい係数を掛けて合計する。求めたい輝度に関係が少なそうな遠いフォトンの影響を下げようということだ。

円錐フィルタは一次関数的に重み付けを行うもので、計算式は本に記載されているものをそのまま使った。$w_{pc}$ は次の通り。

$$w_{pc}=1-\frac{d_p}{kr}$$

$k$の値により、フィルタの効き方が変わってくるが、式から考えると$k$は小さくても1.0まで、無限大にするとフィルタの効果が消える。例の本には$k$が1.1としてあった。フィルタを導入するにあたり、プログラムの関係する部分を少し変更した。(ちなみに全ソースはこちら)

Tracer.hs

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個のフォトンを抽出した上で、推定に使うかどうかは交点からの距離も条件に追加した。

Tracer.hs

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と書く。

cabalでのインストールは下記のようにすれば良い。

$ cabal install mersenne-random

フォトン放射の際に方向ベクトルを生成するが、その関数で使う乱数ライブラリを取り替えて比較しよう。まずは標準乱数ライブラリ。

Algebra.hs

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乱数。

Algebra.hs

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乱数を使うことにしよう。今回唯一の「成果」かな・・・。

まとめ

今回はよりまともに見える画像を生成するよういろいろ試してみたわけだが、どれもあまりうまくいかなかった。結論としては、

このシーンにとっては放射するフォトン数が少ないこと
大域フォトンマップのみではフォトンが全く届かない領域を扱えない

ということかなと。となれば、俄然（鏡面、拡散）反射に対応してその効果を確かめたい！のだが、レイトレーシングの回も結構続いたので一旦休憩しよう。ちょっと他のネタをやって、また本プログラムの拡張に取り組みたいと思う。