この記事はCalendar for レイトレーシング(レイトレ) | Advent Calendar 2021 - Qiitaの記事として作成され、NEEとMISの実装 (NEE編) - MochiMochi3Dの続きです。

MIS(Multiple Importance Sampling)

　NEEは実用的ではありましたがPathTraceの完全上位互換というわけではなく、光源近くではノイズが強くなることや指向性のあるBSDFに弱いという弱点を持っています。一方でPathTraceは全体的なノイズはNEEより多いのですが、PathTraceは指向性のあるBSDFなどに対しては強いという性質を持ちます。NEEとPathtraceはある意味トレードオフのような関係にあります。

　通常であれば、途中でNEEからPathTraceへとサンプリング戦略を変更することは具体的にどんな時に切り替えるべきかの判定が難しいため、サンプリング戦略の切り替えというのはできないと考えられます。ですが、やはりお互いの弱点をカバーできるような仕組みを考えれば一般的なシーンで通用するレンダリング手法が作れるはずです。

　こうした問題を解決するため、MIS（Multiple Importance Sampling)と呼ばれる手法が作られました。

サンプリング戦略の合成

　ある積分のモンテカルロ積分でいくつかサンプリングする方法があったとします。サンプリングする方法をサンプリング戦略といい、 $i$ 個のサンプリング戦略があるものとして考えます。

$\displaystyle{\int_{\Omega} f(x) d\mu(x)} \tag{1}$

通常であれば１つだけのサンプリング戦略を使って積分をモンテカルロ積分していきます。iサンプリング戦略によるモンテカルロ積分は

$\displaystyle{\frac{1}{N_i} \sum_{j}^{N_i} \frac{f(x_j)}{p_i(x_j)} } \tag{2}$

というように行います。一般のモンテカルロ積分では被積分関数 $f$ が式として与えられており、多くの場合かなり効率的なサンプリング戦略を導くことができるため、１つのサンプリング戦略だけで十分なことが多いです。しかしながら、レンダリング方程式においては被積分関数が今レンダリングしたいシーンに依存し、被積分関数が「実際やるまでわからない」です。

　こうした被積分関数が未知である以上、当然万能で効率的な重点的サンプリングは解析的に求めることは不可能であり、PathTraceやNEEのようにある場合には強いがある場合には弱いといった効率的ではあるが万能ではないサンプリング戦略を使っていくしかありません。

　ですが、PathTraceとNEEは丁度お互いがトレードオフのような関係にあり、もしこの２つのサンプリング戦略でそれぞれサンプリングして出てきた結果を上手いこと合わせることができるのであれば、お互いの長所を合わせることができるのではないかと考えられます。このような考えから生まれたのが複数のサンプリング戦略を合成することでより効率的なサンプリングをするMIS(Multi Importance Sampling)という手法になります。

平均による合成

　複数のサンプリング戦略を合成するといってもどのようにすればいいかはあんまり浮かびません。ぱっと思いつく方法としては２つで別々にサンプリングして、結果を平均するということが考えられます。

$\displaystyle{\frac{1}{N_i} \sum_{j}^{N_i} \frac{f(x_j)}{p_i(x_j)} } \tag{3}$

　実際にやってみますとこのような感じになりました。中央はNEE、右がPT、そしてNEEとPTの結果を平均したものは一番左となっています。

f:id:kinakomoti321:20211223031749p:plain — NEEとPTの平均合成

f:id:kinakomoti321:20211223031136p:plain — NEE 50samplingとPT 100sampling

これを見ると、まさにNEEとPTの特徴が平均化されているような画像です。NEEの弱点である光源付近やハイライト部分が完璧でないもののある程度ノイズは消えています。しかしながら、全体的にノイズが多く、これであればNEE単体の方がまだいいのではとも思えます。しかも、PathTraceとNEEをそれぞれ行っているため、非常に処理時間がかかります。５０サンプリングで行いましたがその処理時間はNEEの３倍近くかかっています。

　そんな時間をかけてもこの程度では正直言ってこれでは使い物になりません。このように単純な平均ではあまり効果的とは言えません。これは片方が収束値に近くとも、もう片方がもし全く異なる値を取った時、その平均は収束値から離れてしまうためです。要するに平均値がノイズがひどい方に引っ張られやすくなってしまいます。これではサンプリング結果の合成はうまくいきません。

　単純な平均ではうまくいかないことはわかりましたが、他の方法であればちゃんとした合成ができる可能性はまだあります。この時合成方法をどうやって探していくかということになりますが、そこで考案されたのがMulti Sample Modelというモデルです。

Multi Sample Model

　複数のサンプリング戦略を扱う場合は単純な平均では不十分でより良い組み合わせ方を考えなくてはなりません。より良い組み合わせ方というのは、「いろいろサンプリングしてその結果の内、収束値に近そうな結果が重視される」ことが必要になります。どうやって収束値に近いかを判定するかはともかく、ウェイトを各結果につけることで「結果の重みづけ」をすることができます。

　ということで、複数のサンプリング戦略で得られたサンプリング結果をウェイトを使用して合成するモデルとしてMulti Sample Modelが考案されました。これは各サンプリング戦略にウェイト関数 $\omega_i$ を取り付け、各結果に重みづけができるようになっています。

$f(x)$ に対してサンプリング戦略が $n$ 個存在しているものとして、multi-sample Modelは

$\displaystyle{ F = \sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i} \frac{\omega_i(X_{i,j}) f(X_{i,j})}{p_i(X_{i,j})}} \tag{2}$

と表されます。各文字について

$n_i$ : i番目サンプリング戦略でのサンプル数,

$X_{i,j}$ : i番目サンプリング戦略でj番目のサンプル（乱数）

$\omega_i$ : i番目サンプリング戦略におけるウェイト関数

$p_i$ : i番目サンプリング戦略におけるPDF

という感じになっています。　これが意味する所はそれぞれのサンプリング戦略で別々に $n_i$ 回サンプリングを行い、その平均を合わせるという形で求めるということである。この際、この合成がうまいこと行くように調整項としてウェイト関数 $\omega_i$ をサンプリングごとに入ってきます。

　Multi Sample Modelはちゃんとウェイト関数が

$\displaystyle{ \sum_{i=1}^{n} \omega_i = 1} \tag{3}$

という条件を満たせば、ちゃんと収束値が積分値に収束します。

バランスヒューリスティック

　こうして複数のサンプリング戦略を合成する手段は式として得られましたが、問題はウェイト関数をいかに決めるかということです。実はMulti Sample Modelにおいて（現実的には）最適解とされるウェイト関数は導出することが可能であり、そのウェイト関数は

$\displaystyle{ \omega_i= \frac{n_i p_i}{\sum_{k=1}^{n}n_k p_k} } \tag{4}$

と表されます。これはバランスヒューリスティックと呼ばれており、Multi Sample Modelでは最も分散が小さくなりますウェイト関数になります。これは意外にも数式的に導出が可能で、私がおまけで書いた記事で導出の流れを見ることができます。

kinakomoti321.hatenablog.com

　これを用いたMulti Sample Modelが複数のサンプリング戦略の結果が一番うまくいく合成方法であることになります。MISはこのバランスヒューリスティックを用いて、複数のサンプリング戦略（PT,NEE)の結果を上手く合成することでより効率的なモンテカルロ積分を行っていきます。

MISの実装

　複数のサンプリング戦略を行うのに必要な理論は整いましたのでここからどうやってMISをやっていくかについて書いていきます。今回の場合はサンプリング戦略としてPathTracerとNEEの２つがあるとします。この時におけるMulti Sampling Modelは

$\displaystyle{ F = \frac{1}{n_{pt}} \sum_{j=1}^{n_{pt}} \frac{\omega_{pt}(X_{pt,j}) f(X_{pt,j})}{p_{pt}(X_{pt,j})} +\frac{1}{n_{nee}} \sum_{j=1}^{n_{nee}} \frac{\omega_{nee}(X_{nee,j}) f(X_{nee,j})}{p_{nee}(X_{nee,j})} } \tag{5}$

と書けます。ウェイト関数は

$\displaystyle{\omega_{pt}(X_{pt,j}) = \frac{n_{pt} p_{pt}(X_{pt,j})}{n_{pt} p_{pt}(X_{pt,j}) + n_{nee} p_{nee}(X_{pt,j})} }$

$\displaystyle{\omega_{nee}(X_{nee,j}) = \frac{n_{nee} p_{nee}(X_{nee,j})}{n_{pt} p_{pt}(X_{nee,j}) + n_{nee} p_{nee}(X_{nee,j})} }$

となります。Multi Sampling Modelではそれぞれのサンプル数は異なってもいいようにサンプル数 $n_{pt}$ と分けていますが同じサンプル数 $N$ で行うこととします。こうすることで総和をまとめると、

$\displaystyle{ F = \frac{1}{N} \sum_{j=1}^{N} (\frac{\omega_{pt}(X_{pt,j}) f(X_{pt,j})}{p_{pt}(X_{pt,j})} +\frac{\omega_{nee}(X_{nee,j}) f(X_{nee,j})}{p_{nee}(X_{nee,j})} ) } \tag{6}$

となります。さらにウェイト関数は

$\displaystyle{\omega_{pt}(X_{pt,j}) = \frac{p_{pt}(X_{pt,j})}{p_{pt}(X_{pt,j}) + p_{nee}(X_{pt,j})} }$

$\displaystyle{\omega_{nee}(X_{nee,j}) = \frac{p_{nee}(X_{nee,j})}{p_{pt}(X_{nee,j}) + p_{nee}(X_{nee,j})} }$

となります。

　(6)式の総和のカッコの意味は１回のサンプリングでPathtraceとNEEをそれぞれ行い、それぞれウェイト関数を掛けたうえで合成した値を寄与とすることを示します。そのため、実はそこまで難しい話ではなく前の記事のNEEの実装に加えてPathTraceで寄与を計算できるようにして、ウェイトさえどうにか計算できれば実装はできるというわけです。

ウェイトの計算

　MISの実装はウェイトさえ計算できれば、ほとんど完成のようなものです。ですがウェイトの計算はいくつか注意すべき点があります。ウェイトの計算の流れをたどりながら説明していきます。

　例としてある時のNEEのパスが作れたとしてNEEのウェイト関数を考えてみます。この時、NEEのPDF[tex:p{nee}]は当然計算できますので問題はありません。しかし、ウェイト関数にはさらにどうゆうわけかPTのPDF[tex:p{pt}]という項が存在します。この項 $p_{pt}$ の意味は「もしPTで方向サンプリングしたとして、このNEEで取った方向ベクトルが得られる確率密度」になります。従って、方向ベクトルからPDFを逆算するような機能が必要になります。（例えばコサイン重点的サンプリングであれば方向の角度から求められます）

f:id:kinakomoti321:20211223020939p:plain — ２つの確率密度

　PT側の確率密度を測定できてもまだ１つ重要な注意点が残っています。それは「方向サンプリングと点サンプリングの確率密度は単位が異なる」ということです。単位が異なる量は当然ながら直接和を取ることができません。なので $p_{pt}$ の単位が点サンプリングに合うように変換を加える必要があります。

　変換はNEEにおける幾何項 $G$ を用いて以下のように行えます。（理由は立体角と面積の変換のため）

$\displaystyle{p_{pt} ' = \frac{G}{|\omega_{xy} \cdot n|} p_{pt}} \tag{7}$

　結果として得られるウェイト関数はこんな感じになります。

$\displaystyle{\omega_{nee} = \frac{p_{nee}}{p_{pt} '+ p_{nee}}}$

　PT側でパスを作る場合もほとんど同様です。[tex:p{pt}]を通常通り得た後、「もし仮にNEEをやった時、その光源の点が選ばれる確率密度[tex:p{nee}]」というのを計算しなくてはなりません。その際も先ほどと同様 $p_{nee}$ は単位が合わないため、変換する必要があります。その変換は先ほどの変換の逆数を掛けてあげればいいわけです。

$\displaystyle{p_{nee} ' = \frac{|\omega_{xy} \cdot n|}{G} p_{nee}} \tag{7}$

　結果として得られるウェイトは

$\displaystyle{\omega_{pt} = \frac{p_{pt}}{p_{nee} '+ p_{pt}}}$

と得られます。

実装例

　ウェイトの計算さえわかればMISはすぐに作ることができます。MIS本体の実装の前に必要なものとしては * 受け取った方向ベクトルのPDFを計算する機能 * 当たった光源点のPDFを計算する機能だけです。前者はBSDFにsamplePDFという関数で実装しており、例としてはLambertはコサイン重点的サンプリングなので

  float samplePDF(const Vec3& wo, const Vec3& wi)const override {
    return std::abs(wi[1]);
  }

としたりしています。後者はMISに直で書いています。

　ではMISの処理の流れを書いていきます

For
レイを飛ばす
if 光源に当たった場合
　　//移動で光源にヒットした場合は寄与を取らずに終了
　　break

// NEEの処理
光源上の点をサンプリングする
衝突点からサンプリングした光源の点へレイを飛ばす
if 光源までに何も衝突しなかった場合
　　output = 寄与の計算
        MISweight = MISのウェイトを計算
　　LTE += output * MISweight

//PathTrace
方向サンプリング
if 光源に当たった場合
         output = 寄与の計算
          MISweight = MISのウェイトを計算
         LET += output * MISweight

//次の衝突点への移動 PathTracerで移動する
方向サンプリング
throughputの計算
次に飛ぶレイの生成

ほとんどNEEと同様の処理であり、別途でPTの処理とMISのウェイトを計算する処理を書くだけでMISが完成します。私の実装を以下に書きます。

    Vec3 integrate(const Ray& ray, const Scene& scene, std::shared_ptr<Sampler> sample) const override {
        float p = 1.0; //Probability Rossian Roulette
        const unsigned int MaxDepth = 100; //Max of RayDepth
        Vec3 LTE = Vec3(0.0); //Result of Radiance
        Vec3 throughput = Vec3(1.0); //ThroughPut of Radiance. Midway Waight
        IntersectInfo info; //information of intersected surface. use as information of Current location
        Ray next_ray = ray; //the ray to shot next

        for (int depth = 0; depth < MaxDepth; depth++) {
            // Russian roulette
            p = std::min(std::max(std::max(throughput[0], throughput[1]), throughput[2]), 1.0f);
            if (sample->sample() > p) break;
            throughput /= p;


            //Ray shot
            IntersectInfo info;
            if (!scene.intersect(next_ray, info)) {
                //Ray into sky
                if (depth == 0) {
                    //first shot
                    LTE = throughput * scene.skyLe(next_ray);
                }
                break;
            }


            if (info.object->hasLight()) {
                //Ray into light
                if (depth == 0) {
                    //first shot
                    LTE += throughput * info.object->Le();
                }
                break;
            }

            //directions
            // wo: 入射方向,wi:反射方向
            Vec3 t, b;
            tangentSpaceBasis(info.normal, t, b);
            Vec3 wo = localToWorld(-next_ray.direction, t, info.normal, b);
            Vec3 wi;
            float pdf;
            Vec3 bsdf;

            //NEE
            {
                wo = -next_ray.direction;
                //光源サンプリング
                IntersectInfo lightInfo;
                Vec3 lightPos = scene.lightPointSampling(sample, lightInfo, pdf);

                Vec3 lightDir = normalize(lightPos - info.position);
                Ray shadowRay(info.position, lightDir);
                lightInfo.distance = norm(lightPos - info.position);

                shadowRay.Max = lightInfo.distance - 0.001f;
                IntersectInfo shadowInfo;

                if (!scene.intersect(shadowRay, shadowInfo)) {
                    float cosine1 = std::abs(dot(info.normal, lightDir));
                    float cosine2 = std::abs(dot(lightInfo.normal, -lightDir));

                    wi = lightDir;

                    Vec3 local_wo = worldtoLocal(wo, t, info.normal, b);
                    Vec3 local_wi = worldtoLocal(wi, t, info.normal, b);

                    bsdf = info.object->evaluateBSDF(local_wo, local_wi);
                    //MISWightの計算
                    float lightPDF = pdf;
                    //幾何項（ヤコビアン）の計算
                    float G = cosine2 / (lightInfo.distance * lightInfo.distance);
                    //PTでのPDF
                    float pathPDF = info.object->BSDFpdf(local_wo, local_wi) * G;
                    float MISweight = lightPDF / (lightPDF + pathPDF);
                    //寄与の追加
                    LTE += throughput * MISweight * (bsdf * G * cosine1 / pdf) * lightInfo.object->Le();
                }
            }

            //Pathtrace
            {
                //方向サンプリング
                float pathPdf;
                wo = worldtoLocal(-next_ray.direction, t, info.normal, b);
                Vec3 bsdf = info.object->sampleBSDF(wo, wi, pathPdf, sample);
                const Vec3 nextDir = localToWorld(wi, t, info.normal, b);

                Ray lightRay(info.position, nextDir);
                IntersectInfo lightInfo;

                if (scene.intersect(lightRay, lightInfo)) {
                    if (lightInfo.object->hasLight()) {
                        //衝突かつその物体がLight
                        float cosine1 = absdot(info.normal, nextDir);
                        float cosine2 = absdot(lightInfo.normal, -nextDir);
                        float dis2 = lightInfo.distance * lightInfo.distance;

                        //幾何項(ヤコビアン)の逆数
                        float invG = dis2 / cosine2;

                        //MISweight
                        //NEEでのPDF
                        float lightPdf = 1.0f / (lightInfo.object->areaShape() * scene.getLightN()) * invG;
                        float MISweight = pathPdf / (pathPdf + lightPdf);

                        //Result
                        LTE += throughput * MISweight * cosine1 * lightInfo.object->Le() * bsdf / pathPdf;
                    }
                }
            }

            //次に進む方向のサンプリング

            wo = worldtoLocal(-next_ray.direction, t, info.normal, b);
            bsdf = info.object->sampleBSDF(wo, wi, pdf, sample);

            const Vec3 next_direction = localToWorld(wi, t, info.normal, b);
            const float cosine = std::abs(dot(info.normal, next_direction));
            throughput *= bsdf * cosine / pdf;

            //次のレイ
            next_ray = Ray(info.position + 0.001f * info.normal, next_direction);
        }

        return LTE;
    }

　この実装のMISでサンプリングした画像は以下のようになりました。

f:id:kinakomoti321:20211223031452p:plain — MIS 50sampling

　画像を比較すると、PTよりノイズは少なく、NEEにあったFireFlyや光源近くのノイズがなくなっており、うまいことお互いの特徴を取って比較的ノイズが少なくなっていることがわかると思います。

終わりに

　MISは意外とNEEの実装をしていればそこまで拡張することなく実装ができるのにもかかわらず、綺麗な画像を出すことができます。実際にBlenderのCyclesも（多分）MISが使われており、NEEに比べ実用性がかなり高い手法となっています。実装してみて試しに500サンプリングぐらいやってみるとリファレンスのような綺麗な画像を作れて、ようやく実用的なレイトレができたんだなぁとうれしかったです。

f:id:kinakomoti321:20211223041604p:plain — MIS 500sampling

　MISはパスのサンプリング方法以外にもどんなモンテカルロ積分でも使える一般的なものであるので他にもいろいろな場所で使われています。例えばBSDFがいくつも合わさった複合BxDFはどのBSDFのサンプリングを使うべきかといった問題があるわけですが、これはすなわちサンプリング戦略がいくつもあることであり、MISを使うことができます（One Sample Modelという形式を使う必要があるはずです）。

　もしよろしければおまけでバランスヒューリスティックの導出についての記事も書きましたので見て頂けると幸いです。

kinakomoti321.hatenablog.com

　Veach板出したかった