Mitsuba 2でinverse renderingしてみました

###はじめに
タイトルにアルファベットが多くて、ちまたで発生しているイキった人みたいになっていますが、決して調子に乗ってるわけではなくて、Mitsuba 2 ( https://www.mitsuba-renderer.org/ )というレンダリングのフレームワークを使って、inverse renderingをしてみたので、その通りのタイトルにしています。

inverse renderingが何かというと、コンピュータグラフィックス分野で、レンダリング結果を目標画像に近づけるために、オブジェクトや光源などの環境を求めるという逆問題があり、それを解く技術をinverse renderingと呼びます。
Mitsuba 2では、このinverse renderingの一機能が組み込まれており、使い勝手が良い技術が実装されているともっぱら噂だったので、興味があり、使ってみた次第です。
ちなみに、例を参考に、コードを書いているだけなので、この記事では、Mitsuba 2を使って、Mitsuba 2に用意された機能を実験しただけの内容です。
ただ、例を動かして喜んでいるだけでは芸がないので、少しだけ発展的なことをやってみます。

このページを読み終えて、windows 10でNVIDIAのGPUが刺さっているPCを持っていれば、以下のような結果を作れるようになります。
まず、左の図のような見た目になる光源環境、オブジェクトのマテリアルが設定されたシーンがあるとして、光源環境は変えないで、マテリアルの反射だけを変えて、右の図のような見た目にしたいと思ったとします。
![asampleout0_000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/97cb53c2-5dd9-4344-8879-3dd5bb9c4afe.png)![a0_sample.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/8e3d1686-1813-14c7-f160-bcb7bbfbe5a0.png)
inverse renderingができれば、マテリアルの設定を自動で修正してくれて、制約の中でも、近い見た目を作ることができます。
一例が下の図です。これは初期状態から色の調整をして、近づけていく様子をアニメーションにしてみました。
![onesample.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/e9c151b8-3f71-3d4a-cae4-d653e105afde.gif)
無論変更できないような部分の反射率は変わらないので、全体がそういう色になるとかではありません。
変えられるところのパラメタを自動で修正するというものです。

このページではインストールから、inverse renderingで結果を出すまでの流れを説明します。

####注意
説明すると言ったそばからですが、ところどころに、「本家サイトを読んでくれぇ」と書いてあります。
本家の解説がわかりやすいので、本当に読めばわかることばかりです。
このページでは、わかりにくいところ、若干はまったところを重点的に書いていきます。

あと、コーディングの簡単さについて軽く触れておきます。
私は普通の人なのでプログラミングはあまり得意ではないですが、一応思ったように動かせる程度のことまではできるくらいに簡単です。
Mitsuba2でinverse renderingするためには、Pythonプログラミングが必要となりますが、まともにPythonを書いた経験がない私でもできる範囲のことだけで、inverse renderingできるようになります。
これまでレンダラーを書いて最適化してということをしていた身としては、かなり衝撃的なことです。

また、先に文体や構成について言い訳します。
私は、普段はアル中のごとく、飲み屋のブログくらいしか読まないため、このような様態の技術記事の書き方を理解していません。今回は、周囲の記事を読み、なんとなく皆と同じ様相を持つ文体で記すことに努めました。
また、Markdownの記述も不慣れなため、構造的に読みにくいところがあってもご容赦いただきたいです。

####私の実験環境
読み終えたあとに「私のGPUじゃ合わないじゃん!」とかならないために、先に私の実験PCの構成を書きます。重要なのは、GPUだけですが、一応OS、CPU、RAMも並べておきます。
- Windows 10 64bit
- Intel Core i7 
- NVIDIA GeForce GTX 1080 Ti
- 64GB RAM

GPUは、NVIDIA製で、CUDAが使えないとダメです。お持ちでない方は買ってくるしかないみたいです。
試験的に別コンピュータで試しましたが、GPUが刺さっていない場合、プログラミング環境を構築できませんでした。

ソフト側の環境は、以下です。
- Visual studio 2019
- Miniconda 3
普通ですね。

###少しだけMitsuba 2のこと
コンピュータグラフィックスの学術界隈で、2010年にWenzelさんがMitsubaという物理則に従ったレンダリングのフレームワークを開発しました。
Mitsubaは正確なレンダリング結果を出力するので、精度に厳しい界隈の研究者らは、Mitsubaを使ってレンダリングするようになってきました。
これによってWenzelさんはすごく有名になったのですが、それだけでは止まらず、さらにいくつかの機能を追加したいと思っていたようです（去年本人がどっかでしゃべってた）。
そして、やりたかったことを含めたMitsuba 2が2019年に開発され、SIGGRAPH Asia 2019で発表されました。
その機能のひとつがinverse renderingのひとつの手法であるdifferentiable renderingでした。

ここでは、このdifferentiable renderingを使ってみます。

###インストール
以下のページを開いて、読んだ通りにすればインストールできました。
https://mitsuba2.readthedocs.io/en/latest/

簡単に流れを説明します。
最初のGetting startedと次のCloning the repositoryは、読んだ通りです。
やってください。

Choosing variantsからは注意深くいきましょう。
必ず、Configuring mitsuba.confの項を読んでください。
"enabled": で、gpu_が付いたものを選択しないと、目的のdifferentialbeできないので後悔することになります。
説明文では、gpu_autodiff_rgb かgpu_autodiff_spectralを入れたほうがいいよと書いてありますが、それだけでも後悔するので、以下をそのままコピーしたほうがいいです。

>
"scalar_rgb",
"scalar_spectral",
"gpu_autodiff_rgb",
"gpu_autodiff_spectral",
"gpu_rgb",
"packet_rgb",
"packet_spectral"

次にCompiling the systemでは、各OSの説明の前に、もうGPUを検索して、GPU variantsに飛んでください。
ここで、CUDA ToolkitとOptiXのインストールを先に済ませておきましょう。
終わったらCMakeして、Mitsuba 2をビルドします。
ビルドは、バッチビルドで、ALL_BUILDを選びましょう。個別にチェックしてもいいですが、面倒です。ここでDegugのチェックなんかつけません。Releaseだけです。
![mitsuba2_build.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/23d279f1-24b4-3afa-0374-9e72d82cf337.png)

ビルドが終わったら、distフォルダの中にdllやexeファイルがあるので、そこを環境変数のPathに入れておきます。
setpath.batを走らせてくれ、と本家サイトでは書いてありますが、管理者権限で実行しても、セキュリティのせいか機能しなかったので、手打ちで入力です。
dllのコピーとかだけだと参照漏れがあったときイラつくので、Pathは頭悪い方法ですが確実です。

インストールはこれで終わりです。
要点は以下です。
-　CUDAインストール
-　OptiXインストール
-　configuration.confでもろもろ追加
-　CMake+ビルド
-　環境変数PathにMitsubaのdistを追加

###Pythonでの使い方
本家サイトにそこそこ書いてあるのですが、初心者の私が、はじめに理解できなかったところを書いていきます。

まずmitsuba 2は、xml形式のシーンファイルを読み込んで、その設定どおりにレンダリングします。

本家に書いていないネタですが、実装されていることとして、カメラのタグであるsensorはxml内でいくつも設定できます。
Pythonを使ったとき、xml内に書かれた順のリストになっています。これを切り替えて、カメラ位置を変えてレンダリングができます。

xmlのなかでは、オブジェクトのメッシュデータまでは記述されていません。
Mitsuba2では、.objとか.plyとか使い勝手の良いデータを読めるようになっているので、それら外部ファイルを読み込んで、メッシュとして扱います。

あとは、これらを読むとできます。
https://mitsuba2.readthedocs.io/en/latest/src/getting_started/file_format.html
https://mitsuba2.readthedocs.io/en/latest/src/plugin_reference/intro.html

###初めてのdifferentiable rendering
とりあえず、球でdifferentiable renderingやってみます。
レンダリングしたら、左図の状態ですが、球の反射率を変更して、中図に近づけて、結果的に右図を作ります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/0913eae0-7f72-d1f4-0678-fc175cd4c1c5.png" width="200" alt="画像A" title="画像A"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/2469e952-61c2-f221-e81d-dd20beada6f5.png" width="200" alt="画像B" title="画像B"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/0235d411-5346-22fb-6614-6a4feba49a12.png" width="200" alt="画像C" title="画像C">

下のようにコードの冒頭で、variantを設定します。ここではdifferentialbe renderingが目的なので、gpu_autodiffを使います。色はrgbとします。ここでエラーが出る人は、インストールからやりなおしです。

```
import numpy
import enoki
import mitsuba
mitsuba.set_variant('gpu_autodiff_rgb')
```

次に、xmlは調べれば誰でもわかるので、いつか別の機会に話すとして、とりあえずなのでコードに球とかライトとかもろもろを組み込んじゃいます。
このページを参考にします。
https://mitsuba2.readthedocs.io/en/latest/src/python_interface/parsing_xml.html

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

from mitsuba.core.xml import load_dict

scene = load_dict({
    "type" : "scene",
    "myintegrator" : {
        "type" : "path",
    },
    "mysensor" : {
        "type" : "perspective",
        "near_clip": 0.10,
        "far_clip": 10.0,
        "to_world" : ScalarTransform4f.look_at(origin=[0, 0, 5],
                                               target=[0, 0, 0],
                                               up=[0, 1, 0]),
        "myfilm" : {
            "type" : "hdrfilm",
            "rfilter" : { "type" : "box"},
            "width" : 512,
            "height" : 512,
        },
        "mysampler" : {
            "type" : "independent",
            "sample_count" : 4,
        },
    },
    "myemitter" : {
        "type" : "point",
        "intensity" : 100.,
        "position" :  [5, 5, 5],
        
    },
    "myshape" : {
        "type" : "sphere",
        "mybsdf" : {
            "type" : "diffuse",
            "reflectance" : {
                "type" : "rgb",
                "value" : [0.5, 0.5, 0.5],
            }
        }
    }
})
```
ここまでで、オブジェクトのもろもろは終わりです。
次にdifferentiableの過程を保存するフォルダを作ります。

```
# output folder
outputfolder = 'output/'
import os
if not os.path.exists(outputfolder):
    os.mkdir(outputfolder)
```

outputというフォルダを作って、その中に画像を保存します。
ここからが、目的のコードです。
まず、何がdifferentiableなパラメタかを確認するために、シーンをチェックします。
そのコードが以下です。

```
# Find differentiable scene parameters
params = traverse(scene)
print(params)
```
これを実行すると、このような文字列が現れます。それぞれシーンの要素で、*が先頭についているものがdifferentialbeです。

>
ParameterMap[
  * PointLight.intensity.value,
    PerspectiveCamera.near_clip,
    PerspectiveCamera.far_clip,
    PerspectiveCamera.focus_distance,
    PerspectiveCamera.shutter_open,
    PerspectiveCamera.shutter_open_time,
    Sphere.to_world,
  * Sphere.bsdf.reflectance.value,
]


ということで、differentialbeのものを選びます。ここでは球の反射率を選択します。
具体的な関数の意味とかはドキュメントに書いてありますので割愛です。

```
# to optimize val
optparam = 'Sphere.bsdf.reflectance.value'

# Discard all parameters except for one we want to differentiate
params.keep([optparam])
print(params)
```

これで反射率だけを最適化するパラメタとして絞ることができました。
次に、レンダリングの画像サイズと目指す画像の読み込みと、最適化対象の反射率の初期化をします。

```
#camera sensor
crop_size = scene.sensors()[0].film().crop_size()

#reference
bitmap_ref = Bitmap('myscene/r.png').convert(Bitmap.PixelFormat.RGB, Struct.Type.Float32, srgb_gamma=False)
image_ref = numpy.array(bitmap_ref).flatten()

# Change parameter 
params[optparam] = [0.5, 0.5, 0.5]
params.update()
```

このあたりで書くのが面倒になってきたので、ざっと行きますが、Adamで、レンダリングの結果と目標の画像を比較して、この差が最小となるように、最適化します。
逐次、結果を先に設定した出力フォルダに保存していきます。

```
# Construct an Adam optimizer that will adjust the parameters 'params'
opt = Adam(params, lr=.02)
time_a = time.time()
iterations = 100
for it in range(iterations):
    # Perform a differentiable rendering of the scene
    image = render(scene, optimizer=opt, unbiased=True, spp=3)
    write_bitmap(outputfolder + 'out%03i.png' % it, image, crop_size)
    ob_val_t = enoki.hsum(enoki.sqr(image - image_ref)) / len(image)
    # Back-propagate errors to input parameters
    enoki.backward(ob_val_t)
    # Optimizer: take a gradient step
    opt.step()
    print('Iteration %03i, %s' % (it, params[optparam]), end='\r')
```

100回ループしたら、最後に、時間を計測して終了です。

```    
time_b = time.time()
print('%f ms per iteration' % (((time_b - time_a) * 1000) / iterations))
```

結果は、以下のような様子です。gifアニメなので画質があれですが、実際はもっと階調が豊かです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/e297009b-bae5-6637-16bf-50fe50268a89.gif" width="200">
どんどん目標画像に向かうのがわかります。

これ、何が起きているんだ？と思うかと思いますが、連続的なパラメタ指定できるものを何でも微分できるようなフレームワークである、autodiffを使って、勾配方向を推定し、勾配法でコストが最小となるようにパラメタを寄せていっています。

###今回の実験設計
さて、上のテストでdifferentiable renderingの使い方が分かったので、次は、単一視点のカメラではなく、複数の視点を与えてみます。
それで、視点ごとに色が変わるようなオブジェクトを作ってみます。

コンピュータグラフィックスの研究で、光の吸収率の異なる要素をオブジェクトに良い感じに分布した状態で視点を変えると、同一オブジェクトでも違う見え方になることが明らかになっています。
極論、ホログラムみたいな効果を出すことができるようになります。
こういうものが発展していくと、たとえば、背景に合わせた画像を出力するマントを作れば、纏った者はステルス状態になり究極のボッチを楽しめたり、ひとつの看板であるにもかかわらず見る方向で表示内容が変わり、災害時にどこから見ても適切な逃げ道を提示するような表示ができたりします。
このようなものの基礎となる状態をdifferentiable renderingで作ってみます。

####実験

オブジェクト表面を（誘電体と訳すのでしょうか）dielectricという透明度のあるマテリアルを設定し、光の吸収率の分布をテクスチャで与えるようなシーンを作ります。
また、視点とそれに対応する画像も与えます。
テクスチャを最適化のパラメタとして、狙いの画像に近づくように最適化することで、視点ごとに異なる模様が見えるようになるはずです。

stanford bunnyモデルで実験してみます。環境マップはmitsuba公式サイトにあるものを用います。以下が初期状態です。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/a3ca139d-ecd2-4e72-a1a2-b346cddba146.png" width="200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/5a790cbb-1b56-540e-7ba3-5038224fc05a.png" width="200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/87a82439-cfd8-758d-2e7a-881aa9c55c9a.png" width="200">

上のレンダリング結果を下のようなRGBの縞模様を目指して、モデルの透明度を変更してみます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/3d9917ae-680d-8e85-40ad-2b72625428fa.png" width="200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/f38a58c6-c4ca-fa02-ab6b-f4115b0e97ca.png" width="200"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/20c48d6e-047b-bf75-a6a3-89b1722692d1.png" width="200">

データの入力とか設定なんかは、もうこの際、説明するのもだるいので皆様の想像にお任せするとして、最適化のループだけ書いておきます。
ここでは、3画像を目指すので、一回の最適化ループの中で、3回のdifferentiable renderingを実行します。
感覚的な説明になりますが、全パラメタを同時に最適化をせずに、パラメタの部分ごとの最適化を繰り返して収束させます。もしかしたら発散するかもしれませんが、それは仕方ありません。どんなときに発散するかの条件も知りません。

```
iterations = 100
for it in range(iterations):
    ob_val = 0.0
    for ind in range(len(image_ref)):
        # Perform a differentiable rendering of the scene
        image = render(scene, optimizer=opt, unbiased=True, spp=3, sensor_index = ind)
        write_bitmap(outputfolder + 'out%i_%03i.png' % (ind, it), image, crop_size)
        ob_val_t = enoki.hsum(enoki.sqr(image - image_ref[ind])) / len(image)
        ob_val = ob_val+ob_val_t
        # Back-propagate errors to input parameters
        enoki.backward(ob_val_t)
        # Optimizer: take a gradient step
        opt.step()
        
        #constraint
        params[optparam] = numpy.where(params[optparam]>1., 1., params[optparam])
        params[optparam] = numpy.where(params[optparam]<0.01, 0.01, params[optparam])

    write_bitmap(outputfolder + 'texture_%03i.png' % it, params[optparam], (im_res[0], im_res[1]))
    print('Iteration %03i, %s' % (it, ob_val), end='\r')
```

ちなみに、image_refには、目標画像となるRGBの縞模様がそれぞれ入っています。
また、テクスチャとして透明度を操作しているので、テクスチャが[0,1]の間に入ってほしいです。そのため、numpyで1から0の間の値に押し込んでいます。

####結果
一つのオブジェクトに対して、透明度テクスチャをパラメタとして、differentiable renderingしてみた結果をここに示します。

まず最適化の過程を見てみます。
下の図は、それぞれ異なる視点のカメラで単一のオブジェクトを見ています。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/74d4d2cf-85ac-ef26-df14-b5f5b96b4194.gif", width = "200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/822c9735-50f7-92b7-9be1-392fa3c059c2.gif", width="200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/ad7b0cd0-9bc6-a7f8-b831-ff830130cecd.gif", width="200">

最適化が進むについれて、RGBの縞模様が現れてくるのがわかります。
各視点では、以下の見た目になります。止めていればいい感じに見えます。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/f252fd6d-b90f-32c8-bf47-0e467e2548f3.png" width = "200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/7a321c69-cfa5-eb14-762a-c0641fb241d0.png", width = "200"> <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/538a68d5-6414-4de1-9870-6400e0bb26df.png", width = "200">

最後に得られた透明度で、視点を回してみます。狙いの角度の時は、それっぽい絵が出てますが、ほかはいまいちです。一瞬を見逃したらもう色が違うので、わけがわかりません。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888396/7f60c6ec-e38f-3c98-f94e-5d306d886124.gif", width="200">

とりあえず、奥側の面に色が付いたり、屈折で光が届くところに色がついていたりして、指定の視点のみで成立するような結果が得られました。
differentiable renderingすごい。
ということで、ここで終わりです。

###終わりに

この記事では、Mitsuba 2でinverse rendering してみました。とは言っても、例を読んで、リファレンスを読んで修正したくらいなので、特別、新しいことはしていません。
これくらいのことなら、健康な方なら誰でも実装できるだろうが、その辺は遊びだと思って許容してくれるとうれしいです。
ちなみに、最後の視点を移動させる操作は、カメラの座標をエクセルで作って、xmlに貼り付けました。やっぱりPythonわからんです。