最新物理ベースレンダラー Mitsuba2を触ってみる (3) 微分可能レンダリング編

# Mitsuba2とは
アカデミック向けのフリーの物理ベースレンダラー。微分可能レンダリングや偏光レンダリングなどの新機能がある。
他に書いた関連記事は以下です。

 - [最新物理ベースレンダラー Mitsuba2を触ってみる(1) 導入編](https://qiita.com/kurilab/items/34089636fa570d9a7939)
 - [最新物理ベースレンダラー Mitsuba2を触ってみる(2) Pythonから動かす編](https://qiita.com/kurilab/items/2ad1dd3d0676c1cf4c74)
 - [最新物理ベースレンダラー Mitsuba2を触ってみる(4) 偏光レンダリング編](https://qiita.com/kurilab/items/c96ee592dbde0416576e)
 - [BRDFと偏光BRDF(pBRDF)の違い](https://qiita.com/kurilab/items/f9e48045bd1b9c130d0b)

Mitsuba2のインストールまで終わっているという前提です。

# 微分可能(Differentiable)レンダリングとは
そのままですが微分情報を伝搬できるレンダリングのことです。
レンダリングのアルゴリズムを関数$f(x)$、入力シーンを$x$、レンダリング結果を$y$としたとき、下記が成り立ちます。

```math
y = f(x)
```

微分可能なレンダラーでは、この関数$f$を微分して$\frac{dy}{dx}$を求められるので、所望とするレンダリング結果$y$を得るためにどのように$x$を変更すれば良いかが自動的に求まります。
これは従来では困難であった、CGにおける複雑な逆問題、すなわちインバースレンダリングが解けるということです。

![autodiff.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/215987/a84b36a6-8349-bc84-5946-2fa9d37b7009.jpeg)

とは言ってもそれで何が出来るかすぐにイメージするのが難しいと思うので、具体例を用いながら説明します。

## 集光設計
RGBの光を、ガラスのパネルを通して、壁に投影する例を考えます。プロジェクターと壁の間にガラスを挟むイメージです。
![caustic.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/215987/2ea3ef7e-e3e4-1087-bfbd-ade617373747.jpeg)

この時、投影したい画像を$y$としてあらかじめ与えておきます。
すると、その所望の画像（目的画像）$y$を得るような屈折を実現するようなガラスパネルの微細形状(height profile)$x$を勾配法（SGD, Adam等）で自動最適化することができます。
従来ではこのような問題は、光学系の設計→評価→フィードバックというループを回してヒューリスティックに最適化したり、専用の光学ソフトを用いて各種制約を入れて最適化するしかなったのですが、レンダリングという枠組みで手軽に、かつ汎用的に最適化が出来ます。

## 反射率の最適化
オブジェクトの反射率（BSDFパラメータ）も目的画像を満たすように最適化できます。
下図は、一番右の画像を目的画像$y$として、シーンのパラメータである左側の壁の反射率を最適化する様子を示したものです。
![inv_render0.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/215987/10a3bf55-7d2d-5e99-928b-4da21370160c.jpeg)

左の壁は赤色が真値なのですが、適当な値で反射率を初期化しても最終的に真値に収束していることがわかります（画像が粗いのはサンプリング数が少ないためです、毎回サンプリングを十分に行ってしまうと、微分可能レンダリングの場合は非常に計算時間が多くなってしまうため、最適化中はあまりサンプリング数を上げません）。
このような反射率やBSDFパラメータの自動最適化は、従来は相互反射の影響などもあり、非常に困難な問題とされてきました。ですので既存研究では大幅な制約を加えたりして何とか解いてきたのですが、微分可能なレンダリングではそれが正攻法で直接解くことができます。

## 光源の最適化

オブジェクトの反射率やBSDFパラメータなどは、基本オブジェクト単位で一様で、パラメータ数が少ないので単純な問題と捉えることもできます。
Mitsuba2の微分可能レンダリングでは反射率など単一のスカラー値だけでなく、環境光源マップなどの高次元データも自動最適化できます。
下図は光源を適当に初期化しても真値に収束する例です。
![inv_render1.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/215987/2d744ba6-8e25-d6b5-163a-6898addefc8d.jpeg)

上が環境光源マップ$x$、下がレンダリング結果を示しています。最初は一様な真っ白だった光源マップが変化していき、それに伴いレンダリング結果が目的画像$y$に近づいていることがわかります。１つだけ画像例が与えられるだけで、クリエイターが調整しないでも同じような適切な光源を自動設定できるということになります。

他にもボリュームや形状など、原理的には自動最適化できます。ただし形状などはまだ未実装です。
著者によれば[まもなく実装される](https://github.com/mitsuba-renderer/mitsuba2/issues/26)ということです。GitHubをこまめにチェックすると良いと思います。

# 微分可能レンダリングを動かす
公式が、上の例で挙げた、反射率の最適化や光源の最適化をするサンプルコードを用意してくれているのでそれを動かしてみます。

[導入編](https://qiita.com/kurilab/items/34089636fa570d9a7939)でも書いたことですが、微分可能レンダリングをするためには、GPU環境である必要があります。また、`mistuba.conf`で指定するバリアントに微分可能な構成を最低一つ含める必要があります。まずは`gpu_autodiff_rgb`で良いと思います。

以下はWindowsでの手順を書きます。Linuxの場合もほとんど同じだと思います。MacOSはGPUレンダリングが動きません。
[Pythonから動かす編](https://qiita.com/kurilab/items/2ad1dd3d0676c1cf4c74)と同じ要領で、まずはルートディレクトリでmitsubaモジュールへのPATHを通してから、`docs\examples\10_inverse_rendering\`のサンプルコードに移動します。

```
mitsuba2> setpath.bat
mitsuba2> cd docs\examples\10_inverse_rendering\
mitsuba2\docs\examples\10_inverse_rendering>
```

## 反射率の最適化
まずは反射率の最適化をやってみます。動かすためにデータを持ってきます。

```
git clone https://github.com/mitsuba-renderer/mitsuba-data.git
```
以前も述べましたが上記のデータのcboxをディレクトリごと`10_inverse_rendering`配下に置いてください。
下記を実行し、画像が次々に生成されれば成功です。

```
python invert_cbox.py
```

コードの中身を解説します。和文でコメントをつけました。

```python
# 単純な逆レンダリングの例：コーネルボックスの目的画像をレンダリングし、
# シーンパラメータの1つ（左壁の反射率）を微分可能なレンダリングと勾配ベースの最適化を使用し最適化する

import enoki as ek
import mitsuba
mitsuba.set_variant('gpu_autodiff_rgb')

from mitsuba.core import Thread, Color3f
from mitsuba.core.xml import load_file
from mitsuba.python.util import traverse
from mitsuba.python.autodiff import render, write_bitmap, Adam
import time

# コーネルボックスを読み込む
Thread.thread().file_resolver().append('cbox')
scene = load_file('cbox/cbox.xml')

# 微分可能なシーンパラメータを探す
params = traverse(scene)

# 微分したいパラメータを除いて、すべてのパラメータを破棄する
params.keep(['red.reflectance.value'])

# 現在の値を出力し、バックアップを生成する
param_ref = Color3f(params['red.reflectance.value'])
print(param_ref)

# 目的画像をレンダリングする（まだ導関数は使用していない）
image_ref = render(scene, spp=8)
crop_size = scene.sensors()[0].film().crop_size()
write_bitmap('out_ref.png', image_ref, crop_size)

# 左側の壁を白に変更する（初期値）
params['red.reflectance.value'] = [.9, .9, .9]
params.update()

# パラメータの最適化手法としてAdamを選択
opt = Adam(params, lr=.2)

time_a = time.time()

# イタレーション数
iterations = 100

for it in range(iterations):
    # シーンの微分可能レンダリングを実行
    image = render(scene, optimizer=opt, unbiased=True, spp=1)
    write_bitmap('out_%03i.png' % it, image, crop_size)

    # ロスを取る（レンダリング画像と目的画像の間のMSE）
    ob_val = ek.hsum(ek.sqr(image - image_ref)) / len(image)

    # 誤差を入力パラメータに逆伝搬する
    ek.backward(ob_val)

    # 勾配のステップを実行
    opt.step()

    # 反射率についてGround Truthの値と比較した結果を出力する
    err_ref = ek.hsum(ek.sqr(param_ref - params['red.reflectance.value']))
    print('Iteration %03i: error=%g' % (it, err_ref[0]), end='\r')

time_b = time.time()

# イタレーション毎の計算時間の平均を出力
print()
print('%f ms per iteration' % (((time_b - time_a) * 1000) / iterations))
```

## 光源の最適化
次に、光源マップの最適化を行います。
サンプルデータを[こちらからダウンロードして解凍してください](http://mitsuba-renderer.org/scenes/bunny.zip)。
博物館のような環境に囲まれた金属製のBunnyで、メッシュファイルやシーン記述ファイルなどが入っていると思います。

下記を実行すると光源の最適化が始まります。

```
python invert_bunny.py
```

こちらもコードの中身をコメントをつけて解説します。


```python
import enoki as ek
import mitsuba
mitsuba.set_variant('gpu_autodiff_rgb')

from mitsuba.core import Float, Thread
from mitsuba.core.xml import load_file
from mitsuba.python.util import traverse
from mitsuba.python.autodiff import render, write_bitmap, Adam
import time

# シーン(bunny)を読み込む
Thread.thread().file_resolver().append('bunny')
scene = load_file('bunny/bunny.xml')

# 微分可能なシーンパラメータを探す
params = traverse(scene)

# バックアップを作成する
param_res = params['my_envmap.resolution']
param_ref = Float(params['my_envmap.data'])

# 微分したいパラメータを除いて、すべてのパラメータを破棄する
params.keep(['my_envmap.data'])

# 目的画像をレンダリングする（まだ導関数は使用していない）
image_ref = render(scene, spp=16)
crop_size = scene.sensors()[0].film().crop_size()
write_bitmap('out_ref.png', image_ref, crop_size)

# 白で統一された環境照明マップに変更
params['my_envmap.data'] = ek.full(Float, 1.0, len(param_ref))
params.update()

# パラメータの最適化手法としてAdamを選択
opt = Adam(params, lr=.02)

time_a = time.time()

# イタレーション数
iterations = 100

for it in range(iterations):
    # シーンの微分可能レンダリングを実行
    image = render(scene, optimizer=opt, unbiased=True, spp=1)
    write_bitmap('out_%03i.png' % it, image, crop_size)
    write_bitmap('envmap_%03i.png' % it, params['my_envmap.data'],
                 (param_res[1], param_res[0]))

    # ロスを取る（レンダリング画像と目的画像の間のMSE）
    ob_val = ek.hsum(ek.sqr(image - image_ref)) / len(image)

    # 誤差を入力パラメータに逆伝搬する
    ek.backward(ob_val)

    # 勾配のステップを実行
    opt.step()

    # 環境光源マップについてGround Truthの値と比較した結果を出力する
    err_ref = ek.hsum(ek.sqr(param_ref - params['my_envmap.data']))
    print('Iteration %03i: error=%g' % (it, err_ref[0]), end='\r')

time_b = time.time()

# イタレーション毎の計算時間の平均を出力
print()
print('%f ms per iteration' % (((time_b - time_a) * 1000) / iterations))
```
my_envmap.dataが環境光源マップのパラメータで、125KのRGBAビットマップ（１次元化済）です。
その125Kのパラメータを同時最適化することができます。


上記２つのサンプルコードを比較するとわかるのですが、導出したいパラメータを変えることで、ほぼ同じコードで最適化ができるようになっており、非常に一般的に記述できるということがわかります。

# まとめ
Mitsuba2で新たに可能になった微分可能レンダリングについて、その原理とアプリケーション、そしてサンプルコードの実行までをまとめました。従来は逆問題として困難であったインバースレンダリングが、非常に簡潔な記述で解くことが出来るようになる（場合が多くなる）ということで、個人的には革命的な技術だと思っています。

次回は同じく新たに追加された、偏光レンダリングの機能について紹介します。
 - [最新物理ベースレンダラー Mitsuba2を触ってみる(4) 偏光レンダリング編](https://qiita.com/kurilab/items/c96ee592dbde0416576e)


以上です。