レイトレーシング(14): 光沢面にトライ

これまでは物体表面での反射として完全拡散反射と完全鏡面反射のみ扱ってきた。ただ現実世界はほとんどが光沢のある表面（光沢面とする）なので、リアルな画像を生成するには光沢面のサポートが必須である。とくにぼやけたハイライトは3DCGの醍醐味(?)である。今回はそれを実装してみる。

## 0. Progressive Photon Mapping法に特化

その前に・・・

前回記事では、従来のPhoton Mapping法を用いたレイトレーシングプログラムに対し「オプションとして」Progressive Photon Mapping法（以下PPM法）による画質向上の機能を追加した。しかし作ってみると、このPPM法はこれまでのプログラムに比べ、かなり容易に表現力の向上が可能になるとわかった。たとえば「ちゃんとしたアンチエリアシング」や今回取り組む光沢面などである。

逆に、PPMでないと光沢面などの効果をシンプルに実装できないことから、今後はPPMを主体とすることにした。これにより、PPMを使うかどうかといったフラグや場合わけは排除し、PPMに特化した実装に変更する。結果として、実装はシンプルながらPPMの繰り返し回数を増やすほど高品質な画像が生成できるようになると期待する。

### 1. 光沢面の考え方

光沢面は完全拡散面と完全鏡面の中間となる表面状態である。

![図1.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/1ebd9716-3390-6dc5-50e4-86b05ea23ba9.png)

光子が物体表面に当たった時に、完全にランダム（拡散反射）でも完全に一方向（鏡面反射）でもなく、ある程度の確率で鏡面反射方向の周囲に分残して反射する。どの程度の確率でどれほど分散するかは表面の粗さ（$roughness$）で表す。これが0なら完全に平滑で鏡面反射、1ならどの方向（交点から半球上のランダムな方向）へも反射しうるとする（ただ後述の通り、完全な拡散反射は諦めた）。

ここで「ハイライト」について言及しておきたい。昔からあるベタな3DCG画像では、床平面上に球体が置かれ、例えば球の右上ぐらいに明るいハイライトがぼんやりとした白い円で表現されている。このハイライトが画像のリアルさをぐっと引き上げる。昔ながらの手法ではこのハイライトを視線、光源方向などから計算で擬似的に求めていたわけだ。しかし実際のハイライトは「光源の映り込み」である。周りの物体が反射して映り込むのとなんら変わらない。ゆえにこれまでの記事の作例では球の上部に四角いエリアライトがそのまま映り込んでいる（完全鏡面反射だから）。

![図1.2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/b3d75be4-2395-f296-fa8e-6f608025a5d0.jpeg)

ではぼんやりしたハイライトはどうやってできるのかというと、物体表面がざらついているために、表面の微小なところでは、ある部分は光源が映り込み、ある部分は映り込まないといったことが起こっていて、それらを総合すると光源がぼんやり映ることになるからだ。

であれば、擬似的に計算でハイライトを「作る」のではなく、光源の映り込みという単純な処理でハイライトを表現できるようにしたい。PPMを使えばそれがうまくできそうだ。

### 2. 微小平面の法線

光沢面をサポートするために、光子/視線が反射する方向を$ roughness $に応じてランダムに生成してやる必要がある。その方法には大きく下の２種類がある（と思う）。

1. 交点から鏡面反射方向のベクトルを中心としてその周囲に広がった方向を生成する
2. 交点の微小平面法線をランダムに生成し、それを基に反射方向を算出する

今回は2で実装することにした。理由は以下のとおり。

* 1の場合ランダムに方向を生成すると、入射方向寄りの方向が生成できない（わけではないが面倒に思った）(図2.2-a)
* 微小平面法線を使えば、反射ベクトルだけでなく屈折ベクトルも作り出せる（図2.2-b）
* マイクロファセット理論との相性がいいのではないか（と勝手に推測）

![図2.1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/df347c67-6246-7c98-50cf-75a4984dd64f.png)

では実際の作り方だ。

1. まず交点の法線ベクトル$ n $と光子/視線の入射ベクトル$ e $から、$ n $と直行する交点平面上のベクトル$ u $,$ v $を外積を使って求める。（図2.3）

![図2.3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/98e47919-4735-f9a7-7440-e8554fa73241.png)
2. 次にランダムな方向ベクトルを作るため、２つの一様乱数 $\xi_1 [0,1]$ , $\xi_2 [0,1]$ を生成して下式を用いて上記直行ベクトルの長さを計算する($x$ , $y$ , $z$ とする)。
(参考: [Realistic Image Synthesis - BRDFs and Direct Lighting -](https://graphics.cg.uni-saarland.de/courses/ris-2018/slides/09_BRDF_LightSampling.pdf) のp.5)

```math
\begin{align}
x &= \cos{(2\pi\xi_1})\sqrt{(1-\xi_2^{\frac{2}{m+1}})} \\
y &= \xi_2^{\frac{1}{m+1}} \\
z &= \sin{(2\pi\xi_1})\sqrt{(1-\xi_2^{\frac{2}{m+1}})}
\end{align}
```

3. 直行ベクトルと各係数から微小平面の法線ベクトル $n'$ を求める。

これで法線ベクトル $ n' $ を生成できた。ところで粗さ $ roughness $ はどこへ行ったのだろうか？$ roughness $ に応じてランダムに生成される$n'$ の拡がり具合が変わるようにするのではなかったか？上記の長さを計算する式に出てくる $m$ に注目してほしい。この $m$ を$ 0 $から $\infty$ まで変化させることにより生成される方向の拡がり具合を調整できる。$m$ が0なら一様に拡散（完全拡散面）、$ \infty $ なら $ y $ 成分は常に $1$ （完全鏡面）となるので、$ roughness [0,1] $ から $m$ を作り出せればよい。

```math
\begin{align}
roughness &= 0 \rightarrow m = \infty \\
roughness &= 1 \rightarrow m = 0
\end{align}
```

さて問題は、上記のような変換ができ、かつ$roughness$の変化にできるだけリニア=見た目になだらかな変化、となるような関数をどう選ぶかである。いろいろ試行錯誤した結果、今回は下の式で生成することにした。

```math
m = (10^6)^{(1 - \sqrt{roughness})}
```

この式の意図は次の通り。

* $roughness$の変化に対し $m$ の変化が逆（増えると減る）なので、揃えるために$1-roughness$ のような形にしたが、微調整のため $\sqrt{roughness}$ を使った。
* $10^6$としたのは「見た目になだらかな変化」となるような冪指数を試した結果。好み。
* $\infty$ は無理だが、$1-\sqrt{roughness}$が$1$なら非常に大きな数字になってほしいので指数関数を使った。

ただしこの式ではどんな$roughness [0,1]$ でも $m=0$にはできないので、完全拡散反射は表現できないが。

これで微小平面法線（$n'$）を求めることができる。あとは交点法線 $n$ の代わりに$n'$ を使って光線追跡するだけで良い。なお、このようにランダムに法線を選んでしまうとその方向からくる光しかレンダリングに使えない（もっとあらゆる方向からくる光を集める必要があるのに）と思うが、そこでPPMが威力を発揮するわけだ。イテレーション回数を1000とすれば、ランダムに1000種類の方向をトレースした結果を総合できる。パストレーシングのように交点でランダムな多数のレイを生成して追跡、さらに先の交点でランダムな多数のレイを生成して・・・とすると追跡するレイが指数関数的に爆発するが、PPMではイテレーション回数で抑えられるので処理負荷を低くできる（ハズ）。最初のレイ、交点での反射方向のレイを確率で一つ選んで生成しているので、たぶん
統計学的にもよい結果が得られる、と期待しよう。

### 3. 実装と作例

まず物体の材質を表すパラメータに$roughness$から生成される$m$（を含めた冪指数）を`powerGlossy`として追加する。レンダリングのたびにこの計算を行うのを避けるためである。この冪指数を計算している関数が`densityPower`である。

```Surface.hs
densityPower :: Double -> Double
densityPower rough = 1.0 / (10.0 ** pw + 1.0)
  where
    pw = 6.0 * (1.0 - sqrt rough)
```

次に、ランダムな方向ベクトルを生成する`distributedNormal`である。２つの乱数`xi1`, `xi2`を生成したあと、上記式に基づき係数`x`,`y`,`z`を求め、最後にまとめて法線ベクトル`nvec'`にしている。

```Physics.hs
distributedNormal :: Direction3 -> Double -> IO Direction3
distributedNormal nvec pow = do
  xi1 <- MT.randomIO :: IO Double    -- horizontal
  xi2 <- MT.randomIO :: IO Double    -- virtical
  let
    phi = 2.0 * pi * xi2
    uvec0 = normalize $ nvec <*> (Vector3 0.00424 1.0 0.00764)
    uvec = case uvec0 of
      Just v  -> v
      Nothing -> fromJust $ normalize $ nvec <*> (Vector3 1.0 0.00424 0.00764)
    vvec = uvec <*> nvec
    xi1' = xi1 ** pow
    rt = sqrt (1.0 - xi1' * xi1')

    x = cos(phi) * rt
    y = xi1'
    z = sin(phi) * rt

    nvec' = x *> uvec + y *> nvec + z *> vvec    -- n'の生成
  case normalize nvec' of
    Just v  -> return v
    Nothing -> return ex3    
```

あとはこれまで単純に交点法線を使っていたところを置き換えてやれば良い。

```Tracer.hs
traceRay :: Screen -> Bool -> V.Vector Object -> V.Vector Light -> Int
  -> PhotonMap -> Double -> Material -> Ray -> IO Radiance
traceRay !scr !uc !objs !lgts !l !pmap !radius !m0 !r@(_, vvec) 
  | l >= max_trace = return radiance0
  | otherwise     = do
    case (calcIntersection r objs) of
      Nothing            -> return radiance0
      Just (t, p, n, m, io) -> do                 -- n: 交点の法線ベクトル

   : (中略)

        nvec' <- if rough sf == 0.0
          then return n
          else distributedNormal n (powerGlossy sf)
        let
          (rdir, cos1) = specularReflection nvec' vvec    --n'を使って反射ベクトルを求める

   : (後略)
```

光沢面の作例を次に示す(図3.1-a,b,c)。それぞれ黄色いプラスチック、金、黄色みのガラスである。各球は左上から$roughness$ が0.0, 0.1, 0.2, 0.3, 0.4、下段が左から0.5, 0.6, 0.7, 0.8,1.0としてある。徐々にハイライトがぼやけているのがわかるだろう。またガラスでは球を透過して集光模様ができているが、これも$roughness$に応じてぼやけている。（イテレーション回数=300）

![図3.1-a.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/147b5242-254e-fc57-5dac-48e787b141b1.jpeg)
![図3.1-b.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/c591669b-df1c-70e0-ca9c-82d350e0b9a0.jpeg)
![図3.1-c.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/2c36a22a-262b-d1e1-88d5-6024f525d736.jpeg)

他の作例も示そう。

![図3.2-a.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/12741e75-a87b-d0a7-fee3-4b1e44067304.jpeg)
![図3.2-b.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/69542/f622042b-53be-f83e-6fe7-c0dfa0276367.jpeg)

(図3.2-a, b)

### 4. まとめ

今回は光沢面を扱えるよう機能拡張した。作例からもわかるようにぼんやりしたハイライトがそれなりの品質で表現できるようになったと思う。特に金属ではとても本物らしくみえる。

一方、この手のハイライトは鏡面反射BRDFとして研究されており、Cook-Torranceモデルがしばしば使われているそうだ（下式）。([物理ベースレンダリングを柔らかく説明してみる（４）](https://qiita.com/emadurandal/items/76348ad118c36317ec5c))

```math
f_{r,s} = \frac{DGF}{4 \langle n,l\rangle\langle n, v \rangle}
```


式中の$D$ は微小面法線分布関数NDFといい、今回追加したランダムな反射方向ベクトルの生成と非常に深い関係がありそう・・・だが詳しく調べていない。また、$G$で表される幾何減衰項も微小平面の凹凸による光の遮蔽を表現しているそうで、よりリアルな画像生成には不可欠と思われるがこれも今後の勉強ネタである。今後研究の上、プログラムに組み入れよう。