パストレがしたい(パストレ編) - ややプログラム紀行

パストレ編ってなんだよ、って話ですが、昨日「パストレがしたい」という名前の記事を書きながら肝心のパストレまで到達できなかったので、今回の記事で補完しようと思います。
前回の記事です。
dora119.hateblo.jp

前回は直射光だけを考慮するDirectIntegratorを実装しました。これを拡張してパストレーサーであるPathIntegratorを実装しますが、とりあえず比較しやすいようにコーネルボックスを描画するようにしときました。

f:id:dora119:20181223233447p:plain

思ったよりぽくなって若干嬉しいです。

BRDFのサンプリング

パストレーサーを実装するときにBRDFのサンプリングをする必要があるので、先にこちらの話をしておこうと思います。
おさらいをすると、 $L_o(p, \omega_o) = \int_{H^2(n)} f_r(p, \omega_o, \omega_i) L_i(p, \omega_i) |\cos \theta_i| \mathrm{d}\omega_i$ を分散を抑えつつ計算するには一般的にMultiple Importance Sampling(MIS)という手法によりBRDF側でサンプリングした結果とライト側でサンプリングした結果をいいかんじに合成するのが効果的だが、前回の実装だとライト側のサンプリングで十分なのでBRDFのサンプリングはしなくてもよいという話でした。なので今のところBRDF::f関数とLight::sampleLi関数しか実装されていません。
そこで今回はBRDF::sampleF関数を実装しましょう。BRDF::sampleF関数は反射角を受け取ったときに、入射角をサンプリングし、その角度とともに対応する確率密度を返す関数です。

単純に球面から一様分布でサンプリングしても良いですが、Importance Samplingにおいて重要なことは被積分関数に刺さるようなサンプリングをする必要がある、つまり確率密度関数が被積分関数の定数倍に近くなればなるほど良いということだったのでBRDF::sampleF関数もそのようなサンプリングを目指します。
今回はpbrtにおいてデフォルトで採用されている $p(\omega) \propto \cos \theta$ となる確率分布を採用します。というのも、レンダリング方程式をよく見ると最後に $\cos \theta_i$ の項があり、当然 $\cos \theta_i$ の値が1に近いほどモンテカルロ法において結果への寄与が大きいからです。

cosine-weighted sampling

前回の記事ではサンプリングまわりの話は全て回避してきましたが、今回はそうもいかないみたいです。
$p(\omega) = c \cos \theta$ ( $c$ は定数)となるようなサンプリングを行いたいわけですが、結論からいうと円から一様分布でサンプリングをし、それを半単球面に投射したものがその分布になるらしいです。
イメージ的には球面上の微小面が $\times \cos \theta$ で円に投影されることを考えると筋が通っていそうな気もしますが、証明はpbrを参考にしてください。(2D Sampling with Multidimensional Transformations)

これによって、円から一様分布でサンプリングできればいいことになりました。間違えがちですが、これは極座標 $(r, \phi)$ について $r$ と $\phi$ を一様分布でサンプリングすればいいわけではありません。*1
確率密度関数が $p(x, y) = c$ ( $c$ は定数)となり、円全体で積分すると1になってほしいので、単位円の面積 $\pi$ より $p(x, y) = \frac{1}{\pi}$ となります。ヤコビアンを考えると $p(r, \phi) = rp(x,y) = \frac{r}{\pi}$ となるので、
$p(r) = \int_0^{2\pi} \frac{r}{\pi} \mathrm{d}\phi = 2r$
$p(\phi \mid r) = \int_0^1 \frac{r}{\pi} \mathrm{d}r = \frac{1}{2\pi}$
が分かります。後はこれらの情報をもとにサンプリングを行いたいわけですが、ここではThe Inversion Method(逆変換法)というものによって $[0,1$ ]の一様分布から求める確率分布にしたがった確率変数をサンプリングしてみようと思います。

逆変換法のアイディアはとても単純で、確率密度関数ではなく累積分布関数、の逆関数を用います。
確率密度関数 $p(x)$ に対し $P(a \leq X \leq b) = \int_a^b p(x)\mathrm{d}x$ でしたが、累積分布関数は $F(x) = P(X \leq x)$ のことを指します。この関数は当然単調増加し値域は $[0,1$ ]なので、累積分布関数の逆関数 $F^{-1}$ が求められれば $[0,1$ ]の一様乱数を用いて $p(x)$ に基づく確率変数をサンプリングすることができます。

今回の場合に適用すると、 $r$ をサンプリングするには $p(r) = 2r$ の累積分布関数が $r^2$ なのでこの逆関数 $\sqrt{X_r}$ を、 $\phi$ をサンプリングするには $p(\phi \mid r) = \frac{1}{2\pi}$ より累積分布関数が $\frac{\phi}{2\pi}$ 、よって逆関数 $2\pi X_\phi$ を求めればよいことが分かります。(ここで $X_r, X_\phi$ は一様変数)

$r, \phi$ をサンプリングできれば後はそれを $x, y$ 座標に変換し、半単位球面に投射するので $z = \sqrt{1 - x^2 - y^2}$ を計算するだけです。

virtual RGB sampleF(const Vector3f& _wo, Vector3f* _wi, float* _pdf, Sampler& _sampler) const {
    float r = std::sqrt(_sampler.generateCanonical());
    float theta = 2 * Pi * _sampler.generateCanonical();
    Vector<float, 2> p{r * std::cos(theta), r * std::sin(theta)};
    float z = std::sqrt(std::max(0.0f, 1 - p.at(0) * p.at(0) - p.at(1) * p.at(1)));
    (*_wi)[0] = p.at(0);
    (*_wi)[1] = z;
    (*_wi)[2] = p.at(1);
            
    if (_wo.at(1) < 0.0) (*_wi)[1] *= -1;
    *_pdf = pdf(_wo, *_wi);      
    return f(_wo, *_wi);
}
virtual float pdf(const Vector3f& _wo, const Vector3f& _wi) const {
    return (_wo.at(1) * _wi.at(1) > 0.0) ? std::abs(_wi.at(1)) / Pi : 0.0;
}

_sampler.generateCanonical()関数は内部的にはstd::generate_canonicalを実行しているだけです。*2

パストレーシング(理論)

いよいよ本丸ですが、まずは理論的な話から入ろうと思います。
前回の手法では直射光のみを考慮していたため、コーネルボックスの画像を見ても球の影となる部分は真っ黒になっていることが分かると思います。これは本来ならば、光源でないあるオブジェクトの点も光源に照らされて光を反射しているため二次的な光源となりえますが、この間接光の効果を無視していたからです。よりリアルな画像生成のためには、間接光に照らされた地点は再び二次的な光源となることで他の地点を照らすという光の伝播をシミュレートする必要があります。
レンダリング方程式を用いて説明すると次のようになります。*3
$L_o(p, \omega_o) = L_e(p, \omega_o) + \int_{H^2(n)} f_r(p, \omega_o, \omega_i) L_i(p, \omega_i) |\cos \theta| \mathrm{d}\omega_i$
上の $L_i(p, \omega)$ は本来他の地点から放射されたRadianceであるため、地点 $p$ から $\omega$ の方向にレイを飛ばしたときの衝突地点を $t(p, \omega)$ と表すことにしたとき $L_i(p, \omega) = L_o(t(p, \omega), -\omega)$ と書けるはずです。*4さらに $L_o$ を $L$ と書くようにすれば、レンダリング方程式は
$L(p, \omega_o) = L_e(p, \omega_o) + \int_{H^2(n)} f_r(p, \omega_o, \omega_i) L(t(p, \omega_i), -\omega_i) |\cos \theta| \mathrm{d}\omega_i$
となります。この再帰的な式は一つ $L$ を展開することで以下のようになります。
$L(p, \omega_0) = L_e(p, \omega_0) + \int f_r(p, \omega_0, \omega_1) \\ \cdot \left[ L_e(t(p, \omega_1), -\omega_1) + \int f_r(t(p, \omega_1), -\omega_1, \omega_2) L(t(t(p, \omega_1), \omega_2), -\omega_2) |\cos \theta_1| \mathrm{d}\omega_2 \right] \\ \cdot |\cos \theta_0| \mathrm{d}\omega_1$
めちゃくちゃ読みづらくてアレなので $p_0 = p, p_1 = t(p, \omega_1), p_2 = t(t(p, \omega_1), \omega_2)$ としてまとめると
$L(p_0, \omega_0) = L_e(p_0, \omega_0)\\ + \int f_r(p_0, \omega_0, \omega_1) L_e(p_1, -\omega_1) |\cos \theta_0| \mathrm{d} \omega_1 \\ + \int \int f_r(p_0, \omega_0, \omega_1) f_r(p_1, -\omega_1, \omega_2) L(p_2, -\omega_2) |\cos \theta_0| |cos \theta_1| \mathrm{d}\omega_1\mathrm{d}\omega_2$
であると分かります。

同様に無限に展開しつづけると、
$T_n(\omega_0, \cdots, \omega_n) = \prod_{i=0}^{n-1} f_r(p_i, -\omega_i, \omega_{i+1}) |\cos \theta_i|$
$P_n(\omega_0, \cdots, \omega_n) = \int \cdots \int L_e(p_n, -\omega_n) T_n(\omega_0, \omega_1, \cdots, \omega_n) \mathrm{d}\omega_1 \cdots \mathrm{d}\omega_n$
とすればレンダリング方程式は $L(p_0, \omega_0) = \sum_{i=0}^\infty P_n(-\omega_0, \omega_1, \cdots, \omega_n)$ となります。*5 *6
この $T_n(\omega_0, \cdots, \omega_n)$ は $p_n$ から発射されたRadianceがどれくらい $p_0$ に届くかをあらわし、throughputと呼ばれるそうです。

この形に落とし込んだことでようやく各 $P_n$ に関してはモンテカルロ法が適用できるようになりましたが、まだ $L(p_0, \omega_0)$ を求めるには無限和を計算する必要があります。そこで一般的にはRussian Rouletteという手法を用います。

Russian Roulette

基本的なアイディアはモンテカルロ法と変わりませんが、無限和 $\sum_{i=0}^\infty P_n$ を求めたいときに確率 $q$ で計算を切り上げ、代わりに残りの和を0とすることを考えます。ただし、計算したときは $1-q$ で割るものとすると、全体の期待値は素直に総和を計算した場合と一致します。
上の総和をある地点で分割して $F_1$ と $F_2$ の足し算と考えてみると、F_1計算時の期待値は $(1-q)\frac{F_1}{1-q} + q\cdot 0 = F_1$ となり、確実に $F_1$ を計算した場合と一致します。これが全ての項に関して言えるので総和の期待値も変わらないそうです。*7

パストレーシング(実装)

ようやく最後のパストレ実装になります。理論的には各 $P_n$ をモンテカルロ法により求めれば良いわけですが、計算を省くために $T_n$ は毎回前回の結果を利用することを考えます。というのも、モンテカルロ法を行うときに確率変数をサンプリングして $T_n(\omega_0, \cdots, \omega_n)$ を計算することになりますが、 $n$ を増やしていくときに毎回 $T_n$ を一から計算するのはあまりにしんどいので、前回の計算値(betaとします)に新しくサンプリングした確率変数を用いて $T_n \times f_r(p_n, -\omega_n, \omega_{n+1} |\cos \theta_n|$ を新しいbetaにすることでサボろうというわけです。*8
また、愚直に実装すると飛ばしたパスがライトに衝突しない限り $L_e$ が0となってしまうため非常に効率が悪いですが、各 $P_n$ を計算するときにパスの最後の頂点だけライトのみからサンプリングすることで確実にライトの光を伝播させることができます、これをNext Event Estimationというそうです。

以下、クソ実装です。

RGB li(const Ray& _ray, const Scene& _scene) const {
    Sampler sampler;

    RGB l{0.0, 0.0, 0.0};
    RGB beta{1.0, 1.0, 1.0};
    Ray ray = _ray;

    for (int bounces = 0; bounces < 5; bounces++) {
        SurfaceInteraction si;
        float t;
        if (!_scene.intersect(ray, &si, &t)) break;

        auto light = _scene.getLight();
        l += beta * estimateDirect(si, light, _scene);

        Vector3f l_wo = si.worldToLocal(-ray.getDir().normalize()), l_wi;
        float pdf;
        RGB f = si.brdf->sampleF(l_wo, &l_wi, &pdf, sampler);
        if (pdf == 0.0) break;
        beta *= f * std::abs(l_wi.dot(Vector3f{0.0, 1.0, 0.0})) / pdf;

        ray = Ray(si.point, si.localToWorld(l_wi));

        if (beta.max() < 1.0 && bounces > 3) {
            float q = std::max(0.05f, 1.0f-beta.max());
            if (sampler.generateCanonical() < q) break;
            beta /= 1.0 - q;
        }
    }

    return l;
}

基本的な流れは

ray変数に入っているレイを飛ばす
オブジェクトと衝突しなかったら終了
Next Event Estimationによりbeta(throughput項)と前回実装した直射光を計算するestimateDirectをかけあわせたものをlにたす
パスの最後の頂点からBRDF::sampleFにより次に進む方向を決め、betaとrayを更新
パスの長さが4以上となったらRussian Rouletteをはじめる*9

となっています。
beta(throughput)の値が小さいほど全体への寄与率は小さくなり、そのパスはこれ以上伸ばしてもあまり価値がないということを意味しているので、その値に応じてRussian Roulettteで死ぬ確率も高くなるような実装になっています。

f:id:dora119:20181223233828p:plain — 上のコードにより生成された画像

間接光もシミュレートされるようになったので全体的に明るく、球の影も薄くなり壁の色を若干反射しています。
というか思ったより反射していて、ライトを明るくしすぎたからかもしれないです。
パストレーシングのデメリットとして画面が全体的にノイジーになりました。サンプリングを工夫することでノイズを減らすことができますが、今回は最短でパストレをしたかったのでこんくらいのノイズになっています。*10
そして、もう一つのデメリットですがパストレーシングは大量にシミュレートして平均をとる必要があるので、この画像でも1ピクセルにつき512回くらいシミュレートしていたと思います。