読者です 読者をやめる 読者になる 読者になる

Haskellでいってみよう

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

レイトレーシング(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乱数を使うことにしよう。今回唯一の「成果」かな・・・。

まとめ

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

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