積分とサンプリングと検定

レイトレ合宿3!!!アドベントカレンダー第十回目の記事です。


こんにちは、レイ、飛ばしてますか?
レイトレと言えば積分、積分と言えばレイトレです。
ここでは、積分とサンプリング、そしてその検定についてのさわりを書いていきます。

レンダリング方程式

「レイトレは光のシミュレーション」とよく言われます。キャッチーな言い方で、全くレイトレを知らない人になんとなく概略を伝えるという場合では有効だと思いますが、一方で正確性に関してはミスリードを誘うあまり良くない言い方ではないかと個人的には思っています。ふわっと「光のシミュレーション」と呼んでしまうと「レイが光子か何かの実在のものを表していて、実在のものである以上、誰が実装しても本質的には同じようなものになる」という誤解を与えかねないということです。もし本当に「実在のレイ」があるのなら、実装によってその動きが全く変わってしまうのはとてもおかしな話です。

レイトレは 「連続した輝度の推定」 と言い切ってしまった方が誤解が少ないのではと考えています。

推定なので、完璧な答え(真値)が得られる場合の方が少ないです。また推定の方法、つまり実装者の力量によって性能が大きく変わります。まずレイトレで結局解きたいものは次のKajiyaのRendering equationです。

$ L_0( x,\boldsymbol{\omega_0}) = L_e(x,\boldsymbol{\omega_0}) + \int_{ \Omega} f_r( x,\boldsymbol{\omega_i}, \boldsymbol{\omega_o}) L_i(x,\boldsymbol{\omega_i}) (\boldsymbol{\omega_i} \cdot \boldsymbol{n}) d \boldsymbol{\omega_i} $

左辺は特定のレイの輝度を表します。そのレイがカメラから出るレイ(プライマリレイ)であったらそれはイメージセンサーに入る輝度なので、それを画素上で積分すればフィルムに焼きこまれる強度になります。右辺の第一項は自己発光する物体の光学特性を表します。最も重要な最後の項はレイが交差した表面からどの程度輝度が出てくるかを計算しています。この項の積分の中には左辺で出てきた「特定のレイの輝度」が入っているので、この式は再帰的な式(連続で展開する式)であることがわかります。

物体の光学特性は、ユーザーが直接指定することなので即座に判明します。問題は第二項目です。どのようにしてこれの答えを得たらいいのでしょうか?もし積分区間 $\Omega$ (半球上)の被積分関数について完璧に素性が判明していれば、もしかしたら解析的に解ける場合があるかもしれません。実際「再帰は見なかったことにして半球全体から一定の輝度が来ていると仮定」した人工的なものがAmbient occlusionなのですが、これは単純な形状でのみ構成されていれば、解析的に解ける場合があります。しかし再帰は無視せず、単純な形状はなく、均一に輝度が来るなんてありえない、ごくごくありふれたシチュエーションに対しては解析解は見つけることはできなさそうです。

解析解がないのなら、なんとかして出来るだけ良い推定値を得る必要があります。

問題設定

いきなり現実のレンダーで使われているような方法論を書いても、枝葉の手法に振り回されてしまう面白い人になってしまうので、「サンプリングによって推定する」ことに関して、単純化した次のような問題を考えてみます。

設定1. HDR画像環境マップの輝度に対して何らかの関数を掛けた結果(Convolution)の推定値を得たい
設定2. 何らかの関数の正体は不明
設定3. 同じサンプル数で得られる結果の誤差は確率的に少なければ少ないほどよい
設定4. 最初に何らかの事前処理を行ってもよい

少し抽象的なのですが、「何らかの関数」に具体例を当てはめるとスッと理解しやすいかと思います。
$ F(u,v)=1$ を与えれば、「HDR画像の輝度の平均値」。
$ F(u,v)=cos(v\pi)$ を与えれば、「HDR画像を環境マップとして見立てた時に原点を通過する光束」。
$ F(u,v)=cos(\theta)cos(v\pi)$ を与えれば、「HDR画像を環境マップとして見立てた時の照度」。
BSDFと法線と出射方向を与えれば、「HDR環境マップからの寄与」を計算することになります。
後述のサンプリングの性能評価では、最も単純な一つ目の「輝度の平均値の推定」で比べることにします。

同じHDR画像でも、Convolutionする対象が変われば結果も変わりますが、何らかの関数は最後まで不明だったとしましょう。 一体どうすれば効率的にこの計算を行うことができるのでしょうか?

grace-new

方法1. 全テクセルをサンプリングしてしまう

HDR画像の全てのテクセルをサンプルしてしまう最も単純な力技です。真値を得られる唯一の方法でもあります。しかし推定値を得るために全てサンプルするのは明らかに無駄です。さらに今回の問題の対象が最初から たまたま 離散的な場所(テクセル)で値をとる形態なので、「全て」を見るこの方法が成立します。しかし現実問題の多くは対象が連続関数なはずなのでこの方法はうまくいきません。

次のグラフは全テクセルをサンプルし終わるまでの輝度平均値の推定値の推移です。

mean_estimate_method1

収束が非常に鈍いのがわかると思います。真値を得られるといっても、全テクセルをサンプリングし終わるまで真値は得られません。まじめに「サンプリング」する方法を考えましょう。

方法2. 一様サンプリング

次に思いつくであろう方法は、一様ランダムにサンプリングしていく方法です。 これも無限回続ければ真値に収束するはずです。

mean_estimate_method2

序盤の収束具合は、全サンプリングに比べてかなりマシになりましたが、それでも漸近は鈍く安定はなかなかしません。これはHDR画像の輝度の分布に由来します。対象となるHDR画像のテクセルを輝度値でソートした後にプロットしたものが次のグラフです。

loglog_hdr

両対数であることに注意してください。つまり極一部のテクセルに高い輝度が集中し、大半のテクセルは低い輝度の範囲に収まっているということになります。これはつまり何も考えずに一様にサンプリングしても、

  • 高い確率で推定結果に対して影響を与えない低い輝度を引く。
    これはConvolutionの推定結果への寄与も低くなる。
  • 低い確率で推定結果に大きな影響を与える高い輝度を引く。
    Convolutionの推定結果への寄与が大きいが分散も大きくなってしまう。

ということを意味します。わかりやすい単純化した例として確率$p$で$1/p$の寄与が得られ、それ以外は寄与は得られないような状況を想定します。このとき期待値と分散は次のようになります。

$ X_i = \begin{cases} 1/p & \text{probability p} \\ 0 & \text{otherwise} \end{cases} \\ E[X] = (1/p)p + 0(1-p) = 1 \\ V[X] = E[X^2] - E[X]^2 = (1/p^2)p - 1 = 1/p - 1 $

1p

期待値は常に1であるにも関わらず$p$が小さくなればなるほど、つまりHDR画像の例のように少ないテクセルに輝度が集中すればするほど分散が大きく、つまり推定の誤差が大きくなってしまうことがわかります。これらのことから、とにかく輝度が高いところをまずは完全に捉えないといけないということが分かると思います。

方法3. 高い輝度値を選択する

極一部のテクセルにとても高い輝度が存在するなら、その高い輝度部分だけ抜き出してそこだけで全体をサンプリングしたことにしてしまえばよいと考えるかもしれません。

$ X’_i = \begin{cases} X_i & \text{X_i > threshold} \\ 0 & \text{otherwise} \end{cases} \\ E’[X] = u(X’_i) $

実際にこの方法論で行った結果が次のグラフです。

mean_estimate_method3

横軸の単位が変わっていることに注意してください。今までよりも少ないサンプル数で落ち着いているように見えます。しかし問題があります。

  1. 閾値はどう決めればよいのか。毎回ユーザーが決めるのか?
  2. $E’[X]$はどうしても$E[X]$とは違う値になってしまう。biasが乗ってしまう。
  3. 閾値で「これは影響がある」と決めたテクセルだけに仮に絞ったとしても、
    かなりの数のテクセルが残る。毎回その全てをサンプリングするのは現実的ではない。

これではどうしようも使い物にならないということがわかりました。

方法4. Median Cut

先の方法を進化させた方法として、Median Cutがあります。詳細は[1]を見て下さい。簡単に説明すると

  1. エネルギーが同じになるように、かつなるべく領域が細長くならないように2分割する。
  2. 2分割された領域に対して1を再帰的に行う。
  3. 最初に決めた分割数まで来たら分割を止める。
  4. 分割された各々の領域の重心をエネルギーの分布から決め、
    その一点に光源があると思い込んでサンプリングを行う。

mediancut

面積を持つ領域のエネルギーが重心一点に集中しているという仮定を置いた時点でbiasが発生してますが、分割すればするほど面積は小さくなるのでbiasは小さくなっていきます。

mean_estimate_method4

分割数分サンプリングしてしまえば、(biasがあるにせよ)積分し切ったということになるので収束も早くなります。かなりよさそうな方法に見えますが次の欠点があります。

  1. biasが乗ってしまっている
  2. 領域の分割数に融通が利かない(2の累乗縛り)
  3. 追加でいくらサンプリングしても最初に決めた分割数を超えた詳細な情報は得ることはできない

方法5. 重点サンプリング

大本命です。詳細は[2]を参照してください。簡単に説明すると

  1. 行毎に輝度からCDF(累積分布関数)を求めます。
  2. 行全体の輝度の和から行自体のPDF(確率密度関数)を決めます。
  3. その行のPDFから列用のCDFを求めます

CDFを可視化したのが次の画像です。

cdf

サンプリングをするときは次のようにします。

  1. 二つの乱数(u,v)を生成します。
  2. 列用のCDFと乱数uを使って列を選択します。ここで二分探索を行います。
  3. 選択された列のCDFと乱数vを使って行を選択します。ここでも二分探索を行います。
  4. 最後に、乱数の端数を使ってバイリニアサンプリングを行います。

mean_estimate_method5

圧倒的です。欠点らしき欠点はありません。

[0,2] Sequence, Sobol

重点サンプリングに使う乱数をただの一様乱数ではなくもう少し数値積分に特化したものに置き換えるとさらに性能が向上します。$[0,2]$Sequence+Sobolの組み合わせが収束と計算量のバランスが良いようで、Mitsubaとpbrt-v2の低次元Samplerのデフォルト実装(LDSampler)になっています。原理や理屈、実装の詳細については[3]参照してください。

検定

サンプリングの理屈を完全に理解し、完璧に実装できればよいのですが、あまりそういうことにはならないと思います。どうしても数学的、実装的に正しいことしていたかを確認する必要があります。しかし解析的に結果を出せないので、単純に数値を比較することができません。レンダリングした結果の目視や、推定量の推移から漸近する値を目測するというのがそこそこ用いられる方法だと思います。ただ真値が1とわかっており、推定量が0.97だったとしてそれは正しく推定できていると言えるのでしょうか?

統計の力を使いましょう。

状況をまとめると次のようになります。

  • 真値は何らかの方法で得ている。
  • 新しく実装されたサンプラーを使ってサンプリングした推定量と分散とサンプル数がすでにある。
  • サンプル数自体はかなり多い
  • 間違いなく真値に収束しないと考えることは難しいということを確かめたい

最後の一文は統計の方言で、統計では「確実に同じ」と言い切ることはできず「違うと言いきることは難しい」のような言い方しかできないためにこのように書かれています。

この問題は中心極限定理を使えば次のように置き換えることができる

  • 標本(サンプルの集合)が十分なサイズがあるので、この標本自体の分布は正規分布と仮定できる
  • 標本から計算される信頼区間の中に真値が入っていることを確認できれば、「デタラメな推定量を出しているとは言い難い」と言える。

標本の標準偏差は、標準誤差$\sqrt{V[X]/N}$が相当します。分散がわかれば標準偏差もわかり、標準偏差がわかればその標準偏差内がどれだけの確率をカバーするのかもわかります。具体的には、95%のカバーは標準偏差$1.96\sigma$あたり、99%のカバーは標準偏差$2.33\sigma$あたりに相当します。このカバー率のことを信頼水準(この標本が正しい結論を導ける標本である確率)と呼びます。

結局、信頼水準99%の信頼区間は次のように算出できます。

$ E[X] - 2.33 \sqrt{\cfrac{V[X]}{N}} \leq \theta \leq E[X] + 2.33 \sqrt{\cfrac{V[X]}{N}} $

この区間に先に求めた真値が入っていればとりあえず明後日の方向を向いたサンプリングはしていなかったことになります。

このように真値が完全にわかっていると問題をとても簡単な話にできます。真値がぼんやりしていて区間推定しか出来ていない場合は、ウェルチのt検定を用いる必要があり問題が複雑になってしまいます。力技でサンプル数を増やせば真値に近づくので時間があるときにでも、「間違いようのないレベルの単純なサンプラー」を用いて真値を求めておく方が個人的にはいいと思っています。

結び

今回の例は、単なる手法レベルの話ではなかったことに注意してください。輝度の分布がきわめて不均一な場合のサンプリング方法について書きましたが、このような状況はHDR環境マップのサンプリングだけに限らずレンダラーの内部ではごく自然にあらゆるところであらわれるので、ぜひ積極的に分散を減らしていきましょう。

参考

[1] Paul Debevec. A median cut algorithm for light probe sampling. In ACM SIGGRAPH 2006 Courses, page 6. ACM, 2006. http://gl.ict.usc.edu/Data/HighResProbes/
[2] Monte Carlo Rendering with Natural Illumination, 1–6.
[3] M. Pharr and G. Humphreys, Physically based rendering. San Francisco, Calif.: Morgan Kaufmann, 2010. pp. 359-378