レイトレーシング(11): 画質向上策２点

## 今回は寄り道

前回から半年も経ってしまった orz

前回で完全拡散面の反射まで実装した。当然次は鏡面反射、と思われたら申し訳ないところだが、反射（と屈折）は手間がかかるので後回し。今回は寄り道として生成画像の画質改善に取り組む。

## 1. 直接光は古典レイトレーシングで

フォトンマッピングでは、放射するフォトン数に限界があり、多数のフォトンから輝度推定しても"まだら模様"のようになってしまう。また、点光源では本来影の縁はくっきりするものだが、残念ながらぼやけてしまう。

例の本（フォトンマッピング、Herik Wann Jensen、オーム社）では、実践的なアルゴリズムとして「直接光はフォトンマップから推定するのではなく古典レイトレーシングで」と書いてある（ようだ）。確かにランダムなフォトンから推定するより数式できっちり計算できる古典レイトレーシングを用いたほうが画質的には有利だろう。

前回までで、フォトンマッピング法の画像と比較するため古典レイトレーシングも実装してきた。これをそのまま流用すればよい。基本的な考え方は次の通り。

* フォトン追跡では、光源から放射された直後のフォトンはマップに記録せず、次の反射フォトンから記録するようにする（一次フォトンの記録は捨てる）。
* レイトレーシングでは、光源から直接届く光は数式から求め、それ以外の間接光はフォトンマップから推定する。

### 実装

ではまずフォトン追跡の実装から。

```Tracer.hs
tracePhoton :: [Object] -> Int -> Photon -> IO [PhotonCache]
tracePhoton os l (wl, r) = do
  let is = calcIntersection r os
  if is == Nothing
    then return []
    else do
      let (p, n, m) = fromJust is
      i <- russianRoulette wl [reflectance m]
      pcs <- if i > 0
        then reflect p n os l wl
        else return []
      if (useClassicForDirect == False || l > 0) && diffuseness m > 0.0
        then return $ ((wl, initRay p (getDir r)) : pcs)
        else return pcs
```

`tracePhoton`に$l$という第二引数を設けている。これは追跡の深さを表すもので、物体に当たって反射するたびに1ずつ増やして反射方向の追跡をするようにした。つまり$l=0$ならまだ反射していないから「光源からの直接光」である。

```Tracer.hs
  if (useClassicForDirect == False || l > 0) && diffuseness m > 0.0
```

これがその判定をしているところ。

* `useClassicForDirect`は、光源からの直接光を古典レイトレーシングで計算するかどうかのフラグ。`False`はこれまでどおりフォトンマップから推定するという意味。
* `l > 0`は前述の通り一回以上反射したフォトンを意味する。
* `diffuseness m > 0.0`はフォトンが衝突した物体表面が、少しでも拡散反射の要素を持っている=フォトンを記録する必要があるという意味。

これらのいずれにも合致しない場合は、そのフォトンの衝突点ではフォトンマップに記録しない。

次にレイトレーシング部分を見てみよう。

```Tracer.hs
traceRay :: Int -> Double -> KT.KdTree Double PhotonInfo -> [Object]
         -> [Light] -> Ray -> IO Radiance
traceRay 10 _ _ _ _ _ = return radiance0
traceRay l pw pmap objs lgts r
  | is == Nothing = return radiance0
  | otherwise     = return (em + di + ii)
  where
    is = calcIntersection r objs
    (p, n, m) = fromJust is
    em = sr_half *> emittance m
    di = if useClassicForDirect
      then brdf m $ foldl (+) radiance0 $ map (getRadianceFromLight objs p n) lg ts
      else radiance0
    ii = estimateRadiance pw pmap (p, n, m)
```

`where`節の後半で`useClassicForDirect`というフラグをチェックしている。ここの仕組みはこうだ。

* `useClassicForDirect`が真の場合は、物体上の点に注ぐ直接光を古典レイトレーシングで計算する。`getRadianceFromLight`関数がそれ。光源の数だけ繰り返している。偽の場合は、古典では輝度を計算しないので値ゼロ(=`radiance0`)としている。
* その点には間接光も入ってくるので、その輝度は従来通りフォトンマップから推定する。(ifの下の`estimateRadiance`)

なおフォトンマップから求められる輝度は、古典レイトレーシングを使う場合には「直接光のフォトンを記録しない」ようにしてあるので、ちゃんと間接的に注ぐ光だけで輝度推定される。古典を使わない場合は直接光のフォトンも記録されているので従来と同じ結果となる。

### 画像比較

では生成された画像を比べてみる。まずは点光源の場合。

![image1.png](https://qiita-image-store.s3.amazonaws.com/0/69542/0c3821ed-2bb4-52e9-9ef0-db0dd2b917fc.png)

効果は一目瞭然だろう。直接光を古典で計算した場合(図2)は奥の壁、床、球の上面が滲みなく滑らかに描画されている。直接光の当たらない天井や影の部分は間接光が支配的なので致し方ないが。注目は球の影。従来はどうしても縁がぼやけてしまったがそこがくっきりしている。点光源はこうでないと。次に平面光源の場合。

![image2.png](https://qiita-image-store.s3.amazonaws.com/0/69542/21e4bea3-0274-b000-ac33-b436881c4b37.png)

こちらも同様にスムーズだ。最後に平行光源の場合。

![image3.png](https://qiita-image-store.s3.amazonaws.com/0/69542/560181ba-1d37-659b-1f60-edfd80ecfd03.png)

図5は壁面の直接光の当たる境界がぼやけてしまっているが、古典を使うとくっきりしたエッジになっている。本来こうあるべきだ。（こちらの画像の場合天窓からの光以外は全て間接光なので、ほとんどマダラが残ったままだが、ここは我慢）

処理時間も比べてみよう。直接光を古典レイトレーシングに任せることで直接光のフォトンを記録しないと書いた。つまりフォトンマップが小さくなるということだ。結果は次の表のとおり。

|図|画像|直接光|フォトン数|処理時間(MM:SS)|
|---|---|---|--:|--:|
|1|点光源|フォトンマップ|318,164|6:50|
|2|点光源|古典レイトレ|117,852|1:22|
|3|面光源|フォトンマップ|288,976|6:14|
|4|面光源|古典レイトレ|89,420|1:09|
|5|平行光源|フォトンマップ|321,142|23:07|
|6|平行光源|古典レイトレ|120,841|2:01|

どうやらフォトンマップが小さいと輝度推定時にフォトン収集の時間が短くて済むようで、圧倒的に処理時間が短くなることがわかった！特に平行光源の例では1/10に短縮できている。なぜこれほど短縮できるのか（もしくはフォトンマップだけだと時間がかかるのか）はよくわからないが。

## 2. えせアンチエリアシング

もう一つの画質改善ネタであるアンチエリアシングに移ろう。アンチエリアシングの詳しい説明はWebを漁れば山ほど出てくるのでそちらを見ていただくとして、今回実装した仕組みを説明しよう。

### 仕組み

レイトレーシングは画面を格子状に分割し、その分割された小さい四角（ピクセル）の中心に向かって光を逆向きに追跡する（追跡レイとする）。各ピクセルの輝度値がある閾値を超えたら、そのピクセルの内側のいくつかの点を追加で追跡し、その輝度を平均する。

節タイトルで「えせ」としたのは、ちゃんとアンチエリアシングしたければ追加する追跡レイはピクセル内の「ランダム」な場所に対して行うべきだからである。規則的に追跡レイを追加するとどうしてもエリアシングの問題が残るらしい。しかし今回は、小難しいことを考えるのがだるかったので実装の
簡便さを優先した。。。時間があったらちゃんと実装しよう。

まず、すべてのピクセルで追跡レイを追加していたら時間ばかりかかって非効率なので追加すべきかどうかをまず判断する。そのため、とあるピクセル(p)とその周りの8つのピクセルについてそれぞれ差を求め、閾値を超えているか判定する。1点でも閾値を超えていたら追加レイを飛ばす。（下図）

![image4.png](https://qiita-image-store.s3.amazonaws.com/0/69542/bad6df54-f0bb-7d54-760a-a80aa83cce76.png)

なお、このプログラムでは差を求める際に輝度を比較するのではなくRGB値、すなわち「画面に出力する色」に変換してから比較している。なぜなら、結局目に見えるところでの色の違いがある場合に追加レイを飛ばしたいから。最初は輝度で試したが、閾値の設定がシーンの明るさ・光の量で変える必要が
あり設定が難しく、実用的でないと判断した。

次に追加レイを飛ばす場所だが、上記の通り規則的にした。ピクセル内の中心以外の4点だ。単純にピクセルを4分割してそれぞれの中心点に向かって追跡している。

![image5.png](https://qiita-image-store.s3.amazonaws.com/0/69542/c9605706-396e-c4b1-66be-a13324e0bf3d.png)

これら4点と最初に計算した中心点の全5点から得られた輝度を平均するのだが、その方法には単純平均するか重み付けするかがある。結果に微妙な差が出るが大きな問題はないと判断し、今回は単純平均する。ランダムに追加レイを飛ばす場合は中心からの距離に差が出ると思うのでちゃんと重み付けした方がよいだろう。

### 実装

アンチエリアシングを実施するかどうかはメインループの最後のところ。

```Main-rt.hs
  :
  forM_ [0..(V.length pixels - 1)] $ \i -> do
    rgb <- smooth antiAliasing tracer pixels i
    putStrLn $ rgbToString rgb
```

`smooth`関数で処理している。引数の`antiAliasing`はアンチエリアシングを実施するかどうかのフラグ。`smooth`は次の通り。

```Antialias.hs
smooth :: Bool -> (Ray -> IO Radiance) -> V.Vector Rgb -> Int -> IO Rgb
smooth False _ ims i = return (ims V.! i)
smooth True tracer ims i
  | isDifferent i ims imgOffset == False = return (ims V.! i)
  | otherwise = do
    l <- retrace tracer i
    return $ avg ((ims V.! i):l)
```

まずフラグが偽なら何もせず最初に求めた色を返す。真の場合、まずは追跡レイを追加するかどうか判断する(`isDifferent`)。追加が必要なければやはりそのまま色を返す。追加が必要となったら再トレースだ(`retrace`)。そこで得られた色を`avg`で平均して返せばOK。

追加が必要かどうかの判断は次の通り。

```Antialias.hs
isDifferent :: Int -> V.Vector Rgb -> [Int] ->  Bool
isDifferent _ _ [] = False
isDifferent p rs (i:is)
  | p' <  0         = isDifferent p rs is
  | p' >= length rs = isDifferent p rs is
  | df              = True
  | otherwise       = isDifferent p rs is
  where
    p' = p + i
    df = diffRgb (rs V.! p) (rs V.! p')
```

少し説明が必要かもしれない。2番目の引数は、すでに計算済みの画面全部の色値(RGB)の配列(Vector)である。着目しているピクセルは第一引数`p`に配列内の位置として与えられる。`p'`が比較対象となる
周りのピクセル（の位置）は、第3引数で与えられている。これは上の方で次のように定義されている。

```Antialias.hs
imgOffset :: [Int]
imgOffset = [-xres-1, -xres, -xres+1, -1, 1, xres-1, xres, xres+1]
```

画面の横方向画素数`xres`を使い、着目点`p`に対するオフセットを示しているわけだ。このようにして、周りの8個の点と比べ、閾値を超えたら直ちに`True`を返すようにしてある。周りの点が画面をはみ出す場合はもちろん無視だ(条件 `p' < 0`と`p' >= length rs`）。

次に再トレース部分を説明しよう。

```Antialias.hs
retrace :: (Ray -> IO Radiance) -> Int -> IO [Rgb]
retrace tracer p = do
  let p' = (fromIntegral (p `div` xres), fromIntegral (p `mod` xres))
  ls <- mapM tracer $ map generateRay' (map (badd p') blur)
  return $ map radianceToRgb ls
```

といっても大したことはない。現在の着目点`p`に対し、1/4ずつずれた4点を求め、それぞれ普通に追跡して結果をリストにして返しているだけだ。ずれた点を計算しやすくするため、`blur`をあらかじめ定義している。

```Antialias.hs
blur :: [(Double, Double)]
blur = [(-0.25, -0.25), (-0.25,  0.25), ( 0.25, -0.25), ( 0.25,  0.25)]
```

### 生成画像の比較

では実際にアンチアリアスを効かせた画像を生成してみよう。従来の画像との比較はつぎのとおり。

![image6.png](https://qiita-image-store.s3.amazonaws.com/0/69542/90e324a9-0ee9-0099-77cd-7dabdb907535.png)
![image7.png](https://qiita-image-store.s3.amazonaws.com/0/69542/8d8a691a-a539-348c-9ccb-653ec96d369e.png)

ぱっと見はよくわからないかもしれないが、特に球の上部や平面光源の縁が気持ちなだらかになったのがわかるだろうか？これがアンチエリアシングの効果だ。ただ、いくつかのピクセルについて追加レイを飛ばすので、当然計算時間が余計にかかる。

|アンチエリアシング|フォトン数|処理時間(MM:SS)|
|---|--:|--:|
|無|89,420|1:09|
|有|89,420|1:52|

画面の複雑さや色合いにもよるので一般にどのぐらい時間が増えるかはなんとも言えないが、この例ではほぼ2倍の計算時間がかかったことになる。この辺は画質とのトレードオフだが、世間のCGソフトは普通アンチエリアシングが効いているからなぁ。

## まとめ

今回は多少なりとも画質を向上する策として、以下の2点の改善を実施した。

* 直接光を古典レイトレーシングで計算する
* アンチエリアシングを実装する

両方とも、簡便な方法を採った割にはそれなりにうまくいったと思う。また、古典レイトレーシングを使うことで大幅な計算時間の短縮ができたことは収穫だった。

次回はとうとう(?)鏡面反射、屈折に対応しよう。うまく行けばいわゆる「集光模様」が得られるはずだ。

---
ここまでのソースは[こちら](https://github.com/eijian/raytracer/tree/version-2.3.0.0)。