この記事はレイトレーシングにおけるGGX分布を使用したMicrofacetのBRDFについてまとめたものです。内容としては

となっています。

※　この記事ではBRDFの計算は接空間上で行うこと、Y-upの左手座標系を使用していることを前提としています。

isotropic（等方性）GGX

Microfacet理論

実際の物体は鏡のような滑らかな鏡面だけではなく、粗い金属面も存在しています。粗い表面では細かい傷といった目に見えないほど小さい凹凸が存在し、それぞれで光が反射することで一方向からの光も様々な方向へと乱反射することで粗い表面の質感というものが表れています。

こうした「粗さ」を持つ鏡面のBSDFを作ろうということで考えられたのがMicrofacetモデルです。MicrofacetモデルではMicrofacetと呼ばれる小さな凹凸が存在するものとして考え、そこに入った光はその小さなMicrofacetで鏡面反射するということを仮定して解析的にBSDFとして計算することができるようにしています。

実際に使用される計算式は以下のようになっています。これはBRDFのものでBTDFはまた別にあります（記事ではBRDFのみ扱います）。

$\displaystyle{ f(\omega_i,\omega_o,n) = \frac{F(\omega_i,\omega_m)G(\omega_i,\omega_o,\omega_m)D(\omega_m)}{4|\omega_o \cdot n| |\omega_i \cdot n|}} \tag{1.1}$

$\omega_i$ : 入射方向ベクトル
$\omega_o$ : 出射方向ベクトル
$\omega_m$ : ハーフベクトル
$n$ : 法線ベクトル

ハーフベクトル $\omega_m$ はジオメトリの法線 $n$ とは別にそこに小さな凹凸があるとして（Microfacet）そのミクロな面における法線を意味します。なので出射方向ベクトルを求める際はこのハーフベクトル $\omega_m$ を法線とした反射を考えます。

f:id:kinakomoti321:20220128043336p:plain — マイクロファセットのあるサーフェイス

　次の項ではBRDFに含まれている関数 $F,G,D$ について述べていきます。また、これら関数は粗さのパラメーターとして0~1の $\alpha$ を持ちます。この $\alpha$ について、単純に0~1を指定するのではなく、roughnessとして0~1の値を設定してそれを二乗した値を使用することが多いです。これは理論的な話ではなく、こうした方が粗さの変化が直感的に近いものになるから使われているだけです（こうしなくても問題はないです）。

$\displaystyle{\alpha = roughness * roughness}$

フレネル項 $F$

これは物体のフレネルを司る項で、角度による反射の度合いを示します。実際の計算ではSchlick近似式と呼ばれる

$\displaystyle{F(\omega_i,\omega_m) = F_0 + (1.0 - F_0)(1.0 - \omega_i \cdot \omega_n)} \tag{1.2}$

という式を使用します。ここでの $F_0$ とはスペキュラF0とも呼ばれているやつでフレネルの最小値（面に対して真正面から見た時のフレネルの値）です。要は金属の色でして、こちら側から決定してあげれば大丈夫です（(0.9,0.9,0.9)みたいに決めてもそこまで問題ないです）。

法線分布関数 $D$

　法線分布関数 $D$ とは与えられた $\omega_m$ に対して、 $\omega_m$ を向いたMicrofacetがどれだけあるかを面積として解析的に表現する関数です。これはつまりMicrofacetがどんな形状化かを法線分布として示す関数となっています。このため、 $D$ は反射の性質を根本的に決める重要な関数となっており、また定義として以下のような性質を持ちます。

$\displaystyle{\int_\Omega D(\omega_m) |\omega_m \cdot n| d\omega_m = 1} \tag{1.3}$

　 $D$ は以上のような定義を満たせばよく、様々な種類があります。有名なものとしてはBeckmann分布といったものがありますが、最も実用上使われるのは題名にもあるようにGGX分布と呼ばれるものです(なぜisotropicというものがついてるのかは後述)。

　GGXの法線分布関数 $D$ の式は以下のようになっています。

$\displaystyle{D(\omega_m) = \frac{1}{\pi \alpha^2 \cos{\theta_m}^4 (1 + \frac{\tan{\theta_m}^2}{\alpha^2})^2 }} \tag{1.4}$

この式の $\theta_m$ はハーフベクトル $\omega_m$ とジオメトリの法線 $n$ がなす角であります。

マスキングシャドウ関数 $G$

Microfacetを考えた際、反射した光が別の凹凸によって遮られる場合というのが考えられます。こうした凹凸の遮蔽を表すのがマスキングシャドウ関数 $G$ となります。当然、マイクロファセットの形状によって遮蔽される形というのは変わりますので $D$ それぞれ専用のGが作られています。

$G$ は一般的にSmithのMasking-shadowing関数が使われており、 $G_1$ という関数で

$\displaystyle{G(\omega_i,\omega_o,\omega_m) = G_1(\omega_i,\omega_m) G_1(\omega_o,\omega_m)} \tag{1.5}$

２つのベクトルそれぞれで計算した積として表現できます。

GGXでは $G_1$ はどのような関数かというと

$\displaystyle{G_1(\omega,\omega_m) = \frac{1}{\Lambda(\omega) + 1}} \tag{1.6}$

という式でして、 $\Lambda$ はラフネスの値 $\alpha$ を用いて、

$\displaystyle{\begin{align} \alpha_G &= \frac{1}{\alpha \tan{\theta}}\\ \Lambda(\omega) &= \frac{-1 + \sqrt{1 + \frac{1}{a_G^2}}}{2} \end{align} } \tag{1.7}$

この際、計算に出てきた $\theta$ は $\omega$ とジオメトリの法線 $n$ がなす角です。

また、 $G$ についてより正確な値が出る計算式としてheight correlated Smith shadowing-maskingというものがあります。これは直接 $\Lambda$ を使って、

$\displaystyle{G(\omega_i,\omega_o,\omega_m) = \frac{1}{\Lambda(\omega_o) + \Lambda(\omega_i) + 1}} \tag{1.8}$

と記述されます。こちらを使っても特に問題はありませんのでどっちを使っても大丈夫です。

重点的サンプリング

　GGXはラフネスの値によって方向性が異なるため、半球一様サンプリングは不適切です。そのため、GGX専用の重点的サンプリングを考えてあげる必要があります。　では、GGXのどの部分に対して重点的サンプリングをするかという問題があります。GGXのBRDFはラフネスが低ければ指向性を持ち、鏡面に近づきます。一方で、ラフネスが高ければ等方性が強まり、ざらざらした質感へと変化します。このようなGGXの特徴の根幹を司る部分というのは法線分布関数 $D$ です。そのため、この $D$ に比例するように重点的サンプリングを考えていく必要があります。

　重点的サンプリングを考える際、pdfとしてその関数は全領域で積分したら１にならなくてはいけません。幸い $D$ は定義から $D |\omega_m \cdot n|$ とすれば、式(1.3)のように積分すると１になる関数です。ですので関数 $D$ |\omega_m \cdot n|をpdfとして扱い重点的サンプリングを考えます。

　また、重要な点として $D$ はハーフベクトル $\omega_m$ の関数のため、 $D$ の重点的サンプリングでは出射方向ベクトル $\omega_o$ ではなく、ハーフベクトル $\omega_m$ をサンプルする形になります。サンプリングは乱数[0,1]を $u,v$ として、偏角 $\theta$ と方位角 $\phi$ を

$\displaystyle{\begin{align} \phi_m &= 2 \pi u\\ \theta_m &= \arctan{(\frac{\alpha \sqrt{v}}{\sqrt{1 - v}})} \end{align} } \tag{1.9}$

という式によって求めます。これによってハーフベクトルの $\omega_m$ の角度が定まったのでxyz座標で

$\omega_m = (\cos{\phi_m} \sin{\theta_m},\cos{\theta_m},\sin{\phi_m}\sin{\theta_m}) \tag{1.10}$

というハーフベクトル $\omega_m$ がサンプリングされます。

f:id:kinakomoti321:20220128234256p:plain — ハーフベクトルのサンプリング

　ですが我々が欲しいのはハーフベクトルではなく出射方向ベクトルです。ここからどうやって出射方向ベクトルを求めるかというと、Microfacet理論ではミクロ的にみれば光はMicrofacet上で鏡面反射していることが前提になっているため、ミクロな法線に相当するハーフベクトルを法線として入射ベクトル $\omega_i$ の反射方向を出射方向ベクトル $\omega_o$ とします。

$\displaystyle{\omega_o = refrect(\omega_i,\omega_m)\tag{1.11}}$

　こうして出射方向ベクトル $\omega_o$ をサンプリングできたのですが、大きな注意点としてpdfにヤコビアンがかかることがあります。

　あくまで(1.9)のサンプリングというのはハーフベクトル $\omega_m$ を作るサンプリングであって、直接的に出射方向ベクトル $\omega_o$ をサンプリングするものではありません。出射方向ベクトル $\omega_o$ はハーフベクトルの反射変換を経て得ています。このため、出射方向ベクトルの確率密度pdfというのは変換を挟む以上 $D|\omega_m \cdot n|$ ではなくなります。

　変換によるpdfの変化というのはハーフベクトルと出射方向ベクトルのヤコビアン $J_m$ によって表され、その値というのは

$\displaystyle{J_m = \frac{1}{4 |\omega_o \cdot \omega_m|}} \tag{1.12}$

となっています。これをハーフベクトルのpdfに掛けてあげれば出射方向ベクトルのpdfというのが得られ、結果としては

$\displaystyle{pdf = D|\omega_m \cdot n| J_m= \frac{D|\omega_m \cdot n|}{4 |\omega_o \cdot \omega_m|}} \tag{1.13}$

というpdfであることに注意してください。　

実装例

　実装の前に $n$ との内積や成す角 $theta$ などが多用されるので関数として以下を定義しておきます

namespace BSDFMath {
  inline float cosTheta(const Vec3& w) { return w[1]; }
  inline float cos2Theta(const Vec3& w) { return w[1] * w[1]; }
  inline float sinTheta(const Vec3& w) { return std::sqrt(std::max(1.0f - cosTheta(w) * cosTheta(w), 0.0f)); }
  inline float tanTheta(const Vec3& w) { return sinTheta(w) / cosTheta(w); }
  inline float tan2Theta(const Vec3& w) { return tanTheta(w) * tanTheta(w); }

  inline float cosPhi(const Vec3& w) {
    float sintheta = sinTheta(w);
    if (sintheta == 0) return 1.0f;
    return std::clamp(w[0] / sintheta, -1.0f, 1.0f);
  }
  inline float sinPhi(const Vec3& w) {
    float sintheta = sinTheta(w);
    if (sintheta == 0) return 0.0f;
    return std::clamp(w[2] / sintheta, -1.0f, 1.0f);
  }

  inline float cosPhi2(const Vec3& w) { return cosPhi(w) * cosPhi(w); }
  inline float sinPhi2(const Vec3& w) { return sinPhi(w) * sinPhi(w); }

}

　これらは接空間で法線が(0,1,0)であることを前提として考えて計算されていることに注意してください。

　では次に $F,G,D$ の実装をしていきます。事前にパラメータとしてF0の値と粗さの $\alpha$ を持っているものとして考えていきます。　まずフレネル $F$ についてはshilick近似で

    Vec3 Fresnel(const float im) const {
        float delta = std::max(1.0f - im, 0.0f);
        return F0 + (Vec3(1.0f) - F0) * std::pow(delta, 5);
    }

と実装できます。

　次にGは以下のように実装できます

    float G1(const Vec3& v)const {
        float delta = alpha * BSDFMath::tanTheta(v);
        float result = 2.0f / (1.0f + std::sqrt(1.0f + delta * delta));
        return result;
    }

    float lambda(const Vec3& v) const {
        float delta = std::max(alpha * BSDFMath::tanTheta(v), 0.0f);
        return std::max((-1.0f + std::sqrt(1.0f + delta * delta)) / 2.0f, 0.0f);
    }
    
    float G2(const Vec3& o, const Vec3& i) const {
        return this->G1(o) * this->G1(i);
    }

height correlated Smith shadowing-maskingの場合は以下のようになります

    float G2(const Vec3& o, const Vec3& i) const {
    \\ height correlated Smith shadowing-masking
        return 1.0f / (1.0f + this->lambda(o) + this->lambda(i));
    }

　Dの実装は以下のようになります。

    float D(const Vec3& m) const {
        const float tan2theta = BSDFMath::tan2Theta(m);
        if (std::isinf(tan2theta)) return 0;
        const float cos4theta = BSDFMath::cos2Theta(m) * BSDFMath::cos2Theta(m);
        const float term = 1.0f + tan2theta / (alpha * alpha);
        return 1.0f / ((PI * alpha * alpha * cos4theta) * term * term);
    }

　また、ハーフベクトルのサンプリングは以下のようになります

    Vec3 samplem(float u, float v) const {
        float theta = std::atan(alpha * std::sqrt(v) / std::sqrt(std::max(1.0f - v, 0.0f)));
        float phi = 2.0f * PI * u;
        return Vec3(std::cos(phi) * std::sin(theta), std::cos(theta),
            std::sin(phi) * std::sin(theta));
    }

　以上の関数が実装をまとめたBRDFの実装は以下のようになります。

Vec3 samplingBSDF(const Vec3& wi, Vec3& wo, float& pdf,
        const std::shared_ptr<Sampler>& sampler) const override {

        Vec3 i = wi;
        Vec3 n = Vec3(0.0, 1.0, 0.0);
        Vec3 m = this->samplem(sampler->sample(), sampler->sample());
        Vec3 o = reflect(wi, m);
        wo = o;
        if (wi[1] < 0.0f) return Vec3(0.0);

        float im = absdot(i, m);
        float in = absdot(i, n);
        float on = absdot(o, n);

        Vec3 F = this->Fresnel(im);
        float G_ = this->G2(i, o);
        float D_ = this->D(m);

        Vec3 brdf = F * G_ * D_ / (4.0f * in * on);

        pdf = D_ * BSDFMath::cosTheta(m) / (4.0f * absdot(m, o));
        return brdf;
    }

出力画像

　 $\alpha$ が小さいほど理想鏡面に近づき、一方で大きいほど粗い質感へと近づきます。 $\alpha = 0.1,0.3,0.5$ でレンダリングした結果は以下のような感じになります。

f:id:kinakomoti321:20220129003512p:plain — a = 0.1,0.3,0.5のレンダリング結果

$\alpha = 0.6,0.8,1.0$ のレンダリング結果は以下のようになります。ラフネスが高いとG項が小さくなるため結構暗めの質感を示すようになります。

f:id:kinakomoti321:20220129004136p:plain — a = 0.6,0.8,1.0のレンダリング結果

anistropic (異方性あり）GGX

　実はisotropic GGXでは表現できない金属の質感というのが存在します。例としては、シンバルや金属のコップなどでよくある直線的な線が目立つような旋盤加工した金属が挙げられます。こうした金属のMicrofacetはレコードの表面のような波状の凹凸として存在しています。

　これが意味することとは「向きによって粗さが違う」ということです。そのため、当然ながら向きによって光の反射具合も変化します。こうした性質を「anistropic（異方性）」と呼びます。isotropic GGXはどんな方向にも一様な粗さを持つものを前提として考えており、こうした理由から旋盤加工した金属の質感は扱うことができないというわけです。

　これに対して、異方性を取り扱えるようにしたのがanistropic (異方性あり) GGXです。anisotropic GGXはu,v方向（接空間における接戦方向それぞれ）でそれぞれ粗さを $\alpha_x,\alpha_y$ と別々に考え、法線分布などを作っています。

　anistropicGGXの見た目に関してはこちらのツイートがわかりやすいです。基本的にハイライトが伸びるような見た目を示します。

異方性反射の基礎知識。
これを知っておくと、異方性に関する操作内容がイメージしやすく、扱いやすくなります。#Blender #b3d pic.twitter.com/haXZpeYuKX
— 灰ならし (@hainarashi) 2021年6月8日

異方性あり法線分布関数、シャドウマスキング関数

　法線分布 $D$ やシャドウマスキング関数 $G$ はisotoropicGGXとそこまで変わらず、 $\alpha_x = \alpha_y$ の時はちゃんとisotropicGGXのものと一致するように作られています。 $D$ については以下のようになります。

$\displaystyle{D(\omega_m) = \frac{1}{\pi \alpha^2 \cos{\theta_m}^4 (1 + \tan{\theta_m}^2 (\frac{\cos{\phi_m}^2}{\alpha_x^2} + \frac{\sin{\phi_m}^2}{\alpha_y^2}) )^2 }} \tag{2.1}$

ここでの $\phi_m$ とは接空間におけるハーフベクトル $\omega_m$ の方位角になります。

　一方でシャドウマスキング関数 $G$ そのものは変える必要はなく、 $\alpha_G$ に含まれる $\alpha$ について、

$\displaystyle{\alpha = \sqrt{\cos{\phi}^2 \alpha_x^2 + \sin{\phi}^2 \alpha_y^2}} \tag{2.2}$

として計算することで計算が可能です。

異方性有り法線分布サンプリング

　 $D$ が変わった以上、重点的サンプリングについてもいくつか変更点があります。 $\theta,\phi$ の式は以下のようになります。

$\displaystyle{\begin{align} \phi_m &= \arctan{(\frac{\alpha_y}{\alpha_x} \tan{2\pi u})}\\ \theta_m &= \arctan{(\sqrt{(\frac{1 - v}{v} (\frac{\cos{\phi_m}^2}{\alpha_x^2} + \frac{\sin{\phi_m}^2}{a_y^2}))^{-1}})} \end{align} } \tag{2.3}$

　ただし、実装上の注意としてはArctanが存在することです。Arctanは多価関数であり、C++では主値として $- \pi / 2 \leq \arctan{x} \leq \pi / 2$ となるようにされています。そのため、乱数 $u$ に応じて出てきた値にオフセットを加える形で数値を合わせる必要があります。実装上の式としては以下のようになります。

$\displaystyle{\begin{align} \phi_m &= \left\{ \begin{array}{} \arctan{(\frac{\alpha_y}{\alpha_x} \tan{2\pi u})} & (u \leq 0.25)\\ \arctan{(\frac{\alpha_y}{\alpha_x} \tan{2\pi u})} + \pi & (0.25 < u < 0.75) \\ \arctan{(\frac{\alpha_y}{\alpha_x} \tan{2\pi u})} + 2\pi & (0.75 \leq u) \end{array} \right. \end{align} } \tag{2.4}$

この点にだけ注意すれば実装は(2.3)を行えば大丈夫です。pdfはisotropicと同様です。

実装例

　異方性有りの場合、粗さのパラメータが $\alpha_x,\alpha_y$ の２つのパラメーターになり、　anisotropic GGXの法線分布関数は

    float D(const Vec3& m) const {
        const float tan2theta = BSDFMath::tan2Theta(m);
        if (std::isinf(tan2theta)) return 0;
        const float cos2phi = BSDFMath::cosPhi2(m);
        const float sin2phi = BSDFMath::sinPhi2(m);
        const float cos4theta = BSDFMath::cos2Theta(m) * BSDFMath::cos2Theta(m);

        float term_alpha = cos2phi / (a_x * a_x) + sin2phi / (a_y * a_y);
        float term = 1.0f + tan2theta * term_alpha;
        return 1.0f / (M_PI * a_x * a_y * cos4theta * term * term);
    }

となります。また、シャドウマスキング関数は $\Lambda$ の $\alpha$ の扱いを変えればいいので、

    float lambda(const Vec3& v) const {
        float alpha = std::sqrt(BSDFMath::cosPhi2(v) * a_x * a_x + BSDFMath::sinPhi2(v) * a_y * a_y);
        float delta = std::max(alpha * BSDFMath::tanTheta(v), 0.0f);
        return std::max((-1.0f + std::sqrt(1.0f + delta * delta)) / 2.0f, 0.0f);
    }

となります。最後にハーフベクトルのサンプリングを以下のように実装すれば異方性GGXの完成です。

    Vec3 samplem(float u, float v) const {
        float tanTerm = std::tan(2.0f * PI * u) * a_y / a_x;
        float phi = std::atan(tanTerm);
        if (u >= 0.75) {
            phi += 2.0f * PI;
        }
        else if (u > 0.25) {
            phi += PI;
        }

        const float cos2phi = std::cos(phi) * std::cos(phi);
        const float sin2phi = std::sin(phi) * std::sin(phi);
        const float term_alpha = cos2phi / (a_x * a_x) + sin2phi / (a_y * a_y);

        const float term = v / ((1.0f - v) * term_alpha);
        const float theta = std::atan(std::sqrt(term));

        return Vec3(std::cos(phi) * std::sin(theta), std::cos(theta),
            std::sin(phi) * std::sin(theta));
    }

出力画像

　anisotropic GGXはハイライトが特定の方向に延びるという特徴があります。a_x < a_yの場合横方向にハイライトが伸び、一方でa_x > a_yの場合は縦方向にハイライトが伸びます。

　実際にレンダリングしてみると以下のような感じになります。

f:id:kinakomoti321:20220128043616p:plain — anisotropicGGX a_x = 0.2 a_y = 0.5

f:id:kinakomoti321:20220128043650p:plain — anisotropicGGX a_x = 0.5 a_y = 0.2

可視法線分布サンプリング

可視法線分布とは

今までやってきた重点的サンプリングは法線分布 $D$ に比例するよう行ってきましたが、BRDFの式(1.1)を見ればわかる通りBRDFの値を決める要素は $D$ だけではなく、マスキングシャドウ関数 $G$ もあります。

法線分布 $D$ からサンプリングをする場合、こうした $G$ を考えないでサンプリングするので折角サンプリングした反射光がGによってかき消されてしまうことがあります。特にラフネスが高い物体の輪郭付近ではG項が小さくなりやすく、この理由から法線分布サンプリングではノイズが比較的表れてしまうという弱点が存在します。

そこで「G項に打ち消されない法線分布」というものを作り、それを元に重点的サンプリングをすればこのようなG由来のノイズが減らすことができます。この「G項に打ち消されない法線分布」は「可視法線分布 $D_\omega$ 」と言い、[Heitz et al. 2014]で可視法線分布サンプリングとして提案されました。

どうやってサンプリングするか

可視法線分布 $D_{\omega_i}$ の定義は以下のようになっています。

$\displaystyle{D_{\omega_i(\omega_m)} = \frac{G_1(\omega_i,\omega_m) |\omega_i \cdot \omega_m| D(\omega_m)}{|\omega_i \cdot n|}} \tag{3.1}$

大まかにみれば法線分布 $D$ とシャドウマスキング関数 $G_1$ がかかっているということでG項を考えた法線分布が単純に作れていることがわかると思います。ちゃんとこれは半球面で全積分するとこれは

$\displaystyle{\int_\Omega D_{\omega_i(\omega_m)} d\omega_m = 1} \tag{3.2}$

１になり、ちゃんと逆関数法でのpdfの条件を満たしています。

では、これを重点的サンプリングするため逆関数を考えていくのですがスロープ表現というものを考えたりと非常に複雑で難しいことを行う必要があります。その方法はtatsyさんの導出から実装までの丁寧な記事がありますのでこちらをご参照ください

tatsy.github.io

この方法は[Heitz et al. 2014b]で提唱されていたものですが、難しくまた一部正確ではなかったらしいです。ですので、新たに[Heitz 2016]でGGXに関してより簡単かつ厳密な可視法線分布サンプリングする方法が提唱されました。こちらはスロープ表現を使用する方法よりもはるかに簡単な手法であり、私はこちらを実装しました。(Beckmannにおける手法も一応あります[Wenzel 2014])

https://hal.archives-ouvertes.fr/hal-01509746/document

サンプリングの手順

[Heitz 2016]の内容を書いていきます。なお論文ではz-upの座標系を用いていましたが、ここの説明ではy-up座標系で行おうと思いますので注意してください。

まず、最初に $\omega_i$ を $\alpha_x$ , $\alpha_y$ で引き延ばしてそれを正規化したものを $V$ とします。これが視点ベクトルに相当します。

$\displaystyle{V=normalize(\alpha_x * {\omega_i}_x,{\omega_i}_y,\alpha_y * {\omega_i}_z)} \tag{3.3}$

この視点ベクトル $V$ から正規直交基底 $T_1,T_2$ を作ります。これの作り方は接空間で作るようなものと同じです。

$\displaystyle{\begin{align} T_1 &= normalize(V \times (0,1,0)\\ T_2 &= normalize(T_1 \times V) \end{align} } \tag{3.4}$

また、パラメーターとして $a$ を以下のように定義します。

$\displaystyle{a = \frac{1}{1+V_y}} \tag{3.5}$

一様乱数 $u,v$ [0,1]を用意し、極座標 $(r,\phi)$ を次のようにサンプリングします。 $\phi$ については $v$ の値によって場合分けがあることに注意してください。

$\displaystyle{\begin{align} r &= \sqrt{u} \\ \phi &= \left\{ \begin{array}{} \frac{v}{a} \pi & (v < a)\\ \frac{v - a}{1 - a} \pi + \pi & otherwise \end{array} \right. \end{align} } \tag{3.6}$

　このようにサンプリングした極座標というのはどこの極座標かというと以下の図で言う青と緑の半円板上での極座標です。青の円盤は $T_2$ に平行な半円盤であり、 $v < a$ の時はこちらからサンプリングします。一方緑の部分は $v > a$ の時にサンプリングする場所となっています。

f:id:kinakomoti321:20220126003839p:plain — [Heitz 2016] P1より引用

　なのでこの極座標よりサンプリングする点 $P$ は

$\displaystyle{ P = \left\{ \begin{array}{} r \cos{\phi} T_1 + r \sin{\phi} T_2 & (v < a)\\ r \cos{\phi} T_1 + r \sin{\phi} V_y T_2 & otherwise \end{array} \right.} \tag{3.7}$

となります。この点 $P$ を視点ベクトル方向へと半球面上の投影し、その点を $N$ とすると $N$ は

$\displaystyle{ N = P + \sqrt{1 - P^2} V }$

と得られます。こうしてできた $N$ について、最初のように異方性でまた伸ばして正規化するとようやくハーフベクトル $\omega_m$ が得られます。

$\displaystyle{ \omega_m = normalize(\alpha_x * N_x,N_y,\alpha_y * N_z) } \tag{3.8}$

以上の手順が可視法線分布で重点的サンプリングする方法になります。あとは同様にハーフベクトルと視点ベクトルから反射ベクトルを求めてあげれば良いです。 pdfについては可視法線から得ているため、（ハーフベクトルのヤコビアンは $J_m$ )

$\displaystyle{ pdf = D_{\omega_i}(\omega_m) J_m }$

となります。

実装例

　実装例は以下のようになります。なおBRDFの計算のマスキングシャドウ関数 $G$ がheight correlatedの方であっても特に問題はありませんのでその辺は気にしなくて大丈夫です。

    Vec3 samplem(const Vec3& V_, float u, float v) const {

        Vec3 V = normalize(Vec3(a_x * V_[0], V_[1], a_y * V_[2]));

        Vec3 n = Vec3(0, 1, 0);
        if (V[1] > 0.99) n = Vec3(1, 0, 0);
        Vec3 T1 = normalize(cross(V, n));
        Vec3 T2 = normalize(cross(T1, V));

        float r = std::sqrt(u);
        float a = 1.0f / (1.0f + V[1]);
        float phi;
        if (a > v) {
            phi = PI * v / a;
        }
        else {
            phi = PI * (v - a) / (1.0f - a) + PI;
        }

        float P1 = r * std::cos(phi);
        float P2 = r * std::sin(phi);
        if (a < v) P2 *= V[1];

        Vec3 N = P1 * T1 + P2 * T2 + std::sqrt(std::max(1.0f - P1 * P1 - P2 * P2, 0.0f)) * V;

        N = normalize(Vec3(a_x * N[0], N[1], a_y * N[2]));

        return N;
    }

法線分布サンプリングとの比較

　法線分布サンプリングが弱いのは $D$ が大きいのにも関わらず $G$ が小さくなる時です。これが起きやすいのは大抵ラフネスが高めで法線に対して入射ベクトルが浅い、つまり球体で言う端っこのあたりです。実際に注意深く見てみると球体の端っこのあたりはノイズが強めです。

f:id:kinakomoti321:20220126020708p:plain — 法線分布サンプリングでのGGX

　そうしたノイズへの対策として可視法線分布サンプリングが考えられたので可視法線分布サンプリングによるGGX球体のノイズを比較してみます。サンプリング回数は共に100でRoughnessの値は0.7として両者によるサンプリング画像は

f:id:kinakomoti321:20220126013851p:plain — 法線分布サンプリングと可視法線分布サンプリング

となりました。端っこの部分だけを注目するとちゃんと可視法線分布サンプリングではちゃんと先ほど述べたノイズが軽減されていることがわかると思います。

ちなみにこのサンプリングによる処理時間は200540msと198181msとほとんど等しいぐらいでした。そのため可視法線分布サンプリングは事実上より実用的なサンプリングとして使うことができると考えられます。Blenderのcyclesでも多分可視法線分布サンプリングを使っているっぽいです（おそらくスロープ空間によるサンプリングを使用しているみたいです）。

デバッグのやり方

BRDFのデバッグは中々これといったエラーを出す方法が少なく難しいですが、GGXについてはいくつかデバッグとして使える方法があるのでいくつか書いておこうと思います

法線分布関数 $D$ のデバッグ

　法線分布関数 $D$ について、定義から $D(\omega_m) |\omega_m \cdot n|$ を半球面で積分すると１になるという性質を持ちます。

そのため、ハーフベクトルをランダムに生成してこの積分をモンテカルロ積分してあげればちゃんと１に近づきます。この性質を利用して $D$ が正しく実装されてるかのテストすることができます。

異方性有りの $D$ のデバッグ

　異方性有り $D$ についてはもちろん $D$ なので上記のデバッグができます。それ以外にも、それぞれのラフネスが一致している $\alpha_x = \alpha_y = \alpha$ 時、 $\alpha$ である等方性 $D$ に完全一致する（式的にも）ため、同一の乱数を与えたりし、内部の数値が等方性の $D$ と一致するかということでもデバッグが可能です。