BRDF(GGX)実装してみた

本記事について

本記事はわたくしが実装したBRDFについて解説するものです。
実装内容の概要紹介であり、BRDF自体についての解説や実装方法の詳細については深く触れていません。そのため、読者の方がBRDFやPBRといったことをある程度知っていることを前提に書いています。

作ったもの

こちらが作成したものです。
ただし、シャドウについてはリアルタイム生成ではなく事前にベイクしたテクスチャを使用しています。

こちらがシャドウをオフにしたものです。こちらの方がBRDFの効果がわかりやすいかもしれません。

今回作成したBRDFについて

物理ベースでのBRDFになります。ディフューズはランバート、スペキュラはマイクロファセット理論に基づく実装となっています。

開発にはDirectX12を使用し、本記事内で記述しているコードはhlslになっています。

ディフューズBRDF

ディフューズ（拡散反射）を表すBRDFについて解説します。

ランバートBRDF

今回、ディフューズBRDFにはランバートモデルを採用しています。
BRDF($f_d$)は以下の通りです。こちらは正規化されたランバート反射となっています。

$$\begin{align} f_d(\bold{k}_d)=\frac{\bold{k}_d}{\pi} \end{align}$$

$\bold{k}_d$	RGBごとの拡散反射率。これは材質によって固有の値になる。

実装

float3 Lambert(
    float3 reflection_rate  // 反射率
)
{
    return reflection_rate / F_PI;  // F_PI は事前に定義された円周率
}

スペキュラBRDF

マイクロファセットモデルにおけるスペキュラBRDFは以下の式で表されます。

$$\begin{align} f_r(\bold{v}, \bold{l}) = \frac{D(\bold{n},\bold{h},\alpha)G(\bold{v},\bold{l},\alpha)F(\bold{v},\bold{h},f_0)}{4(\bold{n}\cdot\bold{v})(\bold{n}\cdot\bold{l})} \end{align}$$

記号	説明
$f_r$	BRDF(Specular)
$D$	法線分布関数
$G$	幾何関数
$F$	フレネル関数
$\alpha$	マイクロサーフェイスの粗さを表すスカラー値
$\bold{v}$	視線ベクトル
$\bold{l}$	光源ベクトル
$\bold{h}$	ハーフベクトル($normalize(\bold{v}+\bold{l})$)
$\bold{n}$	ジオメトリの法線ベクトル
$f_0$	ジオメトリ法線ベクトルと視線ベクトルのなす角が0度のときのフレネル反射率

法線分布関数（GGX）

法線分布関数にはGGXを採用します。

$$\begin{align} D_{GGX}(\bold{n},\bold{h},\alpha)=\frac{\alpha^2}{\pi((\bold{n}\cdot\bold{h})^2(\alpha^2 -1)+1)^2} \end{align}$$

実装は以下

float D_GGX(
    float a,    // alpha, 通常はパラメータ値 roughness の二乗となる
    float NoH   // ジオメトリ法線とハーフベクトルの内積
)
{
    // 変数名は適当なので深く考えないでほしい      
    NoH = saturate(NoH);
    float a2 = a * a;
    float b = 1 + (NoH * NoH) * (a2 - 1);
    return a2 / (b * b * F_PI);
}

roughnessについて

一般的に、GGXでの法線分布を計算する際にラフネス値はユーザー入力値の2乗を使って計算することが多いようです。
これは、計算結果と入力値が線形に対応しているように見せかけるための工夫です。線形（のよう）になっていると見た目の調整が楽になります。

詳しくは↓の28pあたり
Background:PhysicsandMathof Shading

幾何関数（SmithCorrelated_GGX）

幾何減衰関数はSmithモデルを使用します。式としては以下のようになります。

$$G(\bold{v},\bold{l},\alpha)=\frac{2(\bold{n}\cdot\bold{l})}{(\bold{n}\cdot\bold{l})+\sqrt{\alpha^2+(1+\alpha^2)(\bold{n}\cdot\bold{l})^2}} \frac{2(\bold{n}\cdot\bold{v})}{(\bold{n}\cdot\bold{v})+\sqrt{\alpha^2+(1-\alpha^2)(\bold{n}\cdot\bold{v})^2}}$$

左辺のうち、右側がシャドウイング（ライトから見えない）、左側がマスキング（カメラから見えない）になります。

↑の$G$に$\frac{1}{4(\bold{n}\cdot\bold{v})(\bold{n}\cdot\bold{l})}$部分を乗算する。こうすると式を単純にすることができます。

$$V(\bold{v},\bold{l},\alpha)=\frac{G(\bold{v},\bold{l}\alpha)}{4(\bold{n}\cdot\bold{v})(\bold{n}\cdot\bold{l})}$$

これを整理すると

$$\begin{align} V(\bold{v},\bold{l},\alpha)=\frac{0.5}{(\bold{n}\cdot\bold{l})\sqrt{(\bold{n}\cdot\bold{v})^2(1-\alpha^2)+\alpha^2}+(\bold{n}\cdot\bold{v})\sqrt{(\bold{n}\cdot\bold{l})^2(1-\alpha^2)+\alpha^2}} \end{align}$$

実装は以下

float G2_SmithCorrelated_GGX
(
    float NoV,
    float NoL,
    float a
)
{
    NoV = abs(NoV) + 1e-5;
    NoL = saturate(NoL);
    float a2 = a * a;
    float GGXV = NoL * sqrt(NoV * NoV * (1.0 - a2) + a2);
    float GGXL = NoV * sqrt(NoL * NoL * (1.0 - a2) + a2);
    return 0.5 / (GGXV + GGXL);
}

フレネル反射関数

フレネル反射関数にはシュリックの近似を使用します。

$$\begin{align} F(\bold{n},\bold{l})\approx f_0+(1-f_0)(1-(\bold{n}\cdot\bold{l})^+)^5 \end{align}$$

実装は以下

float3 Schlick(
    float3 f0,
    float NoL
)
{
    float a = 1 - saturate(NoL);
    return f0 + (float3(1, 1, 1) - f0) * (a * a * a * a * a);
}

$f_0$の導出

$f_0$の導出はFilamentの実装を参考にしました。
だいたい次のような感じ

• 導体(金属)と誘電体(非金属)を分けて考える
• 導体の場合はベースカラーをそのまま$f_0$値として使用
• 誘電体の場合は反射率(refrectance)から計算

誘電体の場合

今回の実装では、誘電体の$f_0$値を次の式で表しています。

$$0.16\cdot reflectance^2$$

$reflectance$は$[0,1]$の範囲で表される値で、反射率を表します。この値は材質によって変わります。
$0.16$は最大フレネル反射率を表しています。これをグラフにすると以下の通り。

Physically Based Rendering in Filament より

一般的な誘電体材質のフレネル反射率はダイヤモンドの0.16を超えることはありません。そのため、$reflectance=1$のとき値が0.16となるようになっています。

導体の場合

導体（金属）の場合はベースカラーをそのままフレネル反射率として用います。導体の場合は拡散反射光はないので、拡散反射率を保持する必要はありません。

金属が拡散反射しないのは、金属内の自由電子にエネルギーが取られるから、らしいです。詳しくは↓。7ページの「Subsurface Scattering」のところ
Background:PhysicsandMathof Shading

スペキュラBRDF のまとめ

これまでの情報をまとめてスペキュラBRDFを実装します。

float3 BRDF_Specular_GGX
(
    float3 N,           // マクロ法線 
    float3 L,           // ライト
    float3 V,           // 視線
    float3 baseColor,   // ベースカラー
    float roughness,    // ラフネス値
    float3 f0          // f0
)
{
    float3 H = normalize(L + V);
    float NoV = dot(N, V);
    float NoL = dot(N, L);
    float NoH = dot(N, H);
    float LoH = dot(L, H);
    
    float a = roughness * roughness;

    // 法線分布
    float D = D_GGX(a, NoH);
    
    // 幾何減衰
    float G = G2_SmithCorrelated_GGX(NoV, NoL, a);
    
    // フレネル
    float3 F = Schlick(f0, NoL);
    
    // スペキュラBRDF
    return (D * G) * F;
}

導体と誘電体の振り分け

導体と誘電体の振り分けにはmetallicという値を使用します。metallicは$[0,1]$の範囲の値で、0が誘電体、1が導体を表します。
基本的には、非金属であればmetallic=0となり、金属であればmetallic=1となります。metallic=0.5のような材質はあまり存在しません。鉄錆などはmetallicが中途半端な値になるかも。

拡散反射光

導体（金属）ではディフューズ反射光はなく、誘電体のみディフューズ反射光が存在します。そのため、線形補間を使ってディフューズ反射光の強さを調整します。

実装は以下のとおりです。

float3 BRDF_Diffuse_Lambert
(
    float3 baseColor,   // ベースカラー
    float metallic      // メタリック
)
{
    //拡散反射（ランバート）
    return Lambert((1.0 - metallic) * baseColor.rgb);
}

スペキュラ反射光

スペキュラ反射光でmetallicに依存するのはフレネル反射成分です。フレネル反射を計算するときに使用する$f_0$値を計算するときにmetallicを使用します。

$f_0$の値

metallicをつかって導体と誘電体を線形補間します。

実装は以下のとおり

float3 Fresnel0
(
    float3 baseColor,
    float metallic,
    float reflectance
)
{
    // 誘電体(Dielectrics)のf0値
    // 誘電体の場合はマテリアルベースカラーに依存しない（無彩色）
    // また、自然界にある材質だとせいぜい0.16（ダイヤモンド）が最大値である
    // よって、0.16を最大値としてreflectanceをリマッピングする
    float3 f0_D = reflectance * reflectance * 0.16;
    
    // 導体(Conductors)のf0値
    // 導体≒金属の場合はベースカラーがそのままf0値となる
    float3 f0_C = baseColor;
    
    // あとは、metallic値によって線形補間する
    return lerp(f0_D, f0_C, metallic);
}

BRDFを使ったライティング実装

これまでの実装をまとめて、ライティングを実装したものが以下になります。

/*
// 以下はマテリアルの基本情報である。設定は適当に
float3 N;                       // ジオメトリ法線
float3 L;                       // ライト方向ベクトル
float3 V;                       // 視線方向ベクトル
float3 H = normalize(V + L);    // ハーフベクトル
float3 baseColor;               // ベースカラー
float roughness;                // ラフネス
float metallic;                 // メタリック
float reflectance;              // 反射率
*/

// ディフューズ反射成分
float3 Fd = BRDF_Diffuse_Lambert(baseColor, metallic);

// スペキュラ反射成分
float3 f0 = Fresnel0(baseColor, metallic, reflectance);
float3 Fr = BRDF_Specular_GGX(
    N,
    L,
    V,
    baseColor,
    roughness,
    f0
);

// BRDF
float3 BRDF = (Fd + Fr);

以上でBRDFの実装は完了です。

最終的に出力するカラー

BRDFの計算はすでに完了していますが、ここまできたら最終的な出力カラーまで計算してしまいましょう。
以下のコードは直前のコードの続きです。

// ライト。値は適当
float3 lit = float3(1.0, 1.0, 1.0);

// 出力カラー
float3 outputColor = saturate(dot(N, L)) * BRDF * lit;

簡単化のためライトを定数として設定しています。これでもBRDFの確認はできますが、絵的に面白い画面を作ろうと思ったらライトの値も複雑になっていきます。

補足

以下は補足的な内容となっています。
各関数の導出など、調べたことを書いてあります。

Kajiyaのレンダリング方程式

Kajiyaのレンダリング方程式は、Hanes T.Kajiya氏が考案した積分方程式で、物理ベースのレンダリングの計算・実装を行う際によく利用されます。
Real Time Rendering(第4版) 11.1式・11.2式より(ちょっとわかりやすく変えてあります)

$$\begin{align} L_o(\bold{p},\bold{v})=L_e(\bold{p},\bold{v})+\int_{\bold{l}\in\Omega}f(\bold{l},\bold{v})L_i(\bold{p},\bold{l})(\bold{n}\cdot\bold{l})^+ d\bold{l} \end{align}$$

$L_o$	出射光
$L_e$	発光
$f$	BRDF
$L_i$	入射光
$\bold{p}$	サーフェイス上の位置ベクトル
$\bold{v}$	視線ベクトル
$\bold{l}$	ライトベクトル
$\Omega$	半球立体角（積分範囲）

$x^+$は$x$を0にクランプするという意味

BRDFを考える上では発光を表す$L_e$は必要ないため省略することが多いです。

$$\begin{align}L_0(\bold{p},\bold{v})=\int_{\bold{l}\in\Omega}f(\bold{l},\bold{v})L_i(\bold{p},\bold{l})(\bold{n}\cdot\bold{l})^+ d\bold{l} \end{align}$$

正規化ランバートの導出

ランバートモデルは以下のように表されます。

$$L_o=L_i\bold{k}_d (\bold{n}\cdot\bold{l})^+$$

レンダリング方程式を使って、エネルギー保存則を満たすよう変形したものが正規化ランバートモデルとなります。
以下、その導出です。

前提として、今回考える拡散反射は入射した光が半球上に一様に反射する完全拡散面での反射とします。
すると、反射光の比率は入力値に関わらず一定となるため、定数で表すことができます。

$$f_d(\bold{l},\bold{v})=K\bold{k}_d$$

入力値は関係なくなるので、省略します。

$$f_d=K\bold{k}_d$$

$K$はBRDFによる反射率、$\bold{k}_d$は材質毎にことなる反射率を表しています。BRDFを考える上で、材質毎に固有の反射率は考えないため、$\bold{k}_d=(1,1,1)$として省略してしまうことにします。

$$f_d=K$$

これをレンダリング方程式に当てはめます。($f(\bold{l},\bold{v})=f_d=K$となる)なお、拡散反射以外の反射は考えません。

$$\begin{align*} L_0(\bold{p},\bold{v})&=\int_{\bold{l}\in\Omega}K L_i(\bold{p},\bold{l})(\bold{n}\cdot\bold{l})^+ d\bold{l} \\ &=K\int_{\bold{l}\in\Omega}L_i(\bold{p},\bold{l})(\bold{n}\cdot\bold{l})^+ d\bold{l} \end{align*}$$

入射光($L_i(\bold{p},\bold{l})$)も(1,1,1)として省略します。

$$L_0(\bold{p},\bold{v})=K\int_{\bold{l}\in\Omega}(\bold{n}\cdot\bold{l})^+ d\bold{l}$$

$\bold{n}$と$\bold{l}$のなす角を$\theta$とすると$(\bold{n}\cdot\bold{l})=\theta$となるので

$$L_0(\bold{p},\bold{v})=K\int_{\bold{l}\in\Omega}\cos\theta d\bold{l}$$

積分範囲は立体角$\omega$であるため、$d\bold{l}$は次のように表せます。
参考：https://youtu.be/aNoEzONgIYo?t=454

$$d\bold{l}=\sin\theta d\theta d\phi$$

こちらを使って積分範囲を修正します。半球であるため、$0\leq\theta\leq\frac{\pi}{2}$、$0\leq\phi\leq2\pi$となります。
先程の式にあてはめると

$$L_0(\bold{p},\bold{v})=K\int_0^{2\pi}\int_0^{\frac{\pi}{2}}\cos\theta \sin\theta d\theta d\phi$$

$\phi$の部分は定数の積分とみなせるので、$\phi$部分の積分は以下のようになります。

$$\begin{align*} \int_0^{2\pi}d\phi&=[\phi]_0^{2\pi} \\ &=2\pi-0 \\ &=2\pi \end{align*}$$

よって、先程の式は次のようになります。

$$\begin{align*} L_o(\bold{p},\bold{v})&=K\int_0^{2\pi}\int_0^{\frac{\pi}{2}}\cos\theta \sin\theta d\theta d\phi \\ &=2\pi K\int_0^{\frac{\pi}{2}}\cos\theta \sin\theta d\theta \end{align*}$$

二倍角の公式から$\cos\theta \sin\theta$の部分をまとめます。

$$\begin{align*} 2\sin\theta\cos\theta &= \sin2\theta \\ \sin\theta\cos\theta &= \frac{\sin2\theta}{2} \end{align*}$$

よって

$$\begin{align*} 2\pi K\int_0^{\frac{\pi}{2}}\cos\theta \sin\theta d\theta &= 2\pi K \int_0^{\frac{\pi}{2}}\frac{\sin2\theta}{2}d\theta \\ &= \pi K \int_0^{\frac{\pi}{2}}\sin2\theta d\theta \\ &= \pi K \bigg[ \frac{-\cos 2\theta}{2} \bigg]_0^{\frac{\pi}{0}} \\ &= \pi K \bigg(\frac{-\cos \pi}{2} - \bigg(\frac{- \cos 0}{2} \bigg) \bigg) \\ &= \frac{\pi K}{2} (1 + 1) \\ &= \pi K \end{align*}$$

以上をまとめると以下になります。

$$\begin{align*} L_o(\bold{p},\bold{v})&=\int_{\bold{l}\in\Omega}K(\bold{n}\cdot\bold{l})^+ d\bold{l} \\ &= \pi K \end{align*}$$

これは、反射率が方向によらず一定($K$)である場合、反射光の強さは$\pi K$になるということを表しています。

入射光$L_i$が$1$であるという前提であったため、出射光$L_o$は$1$を超えてはいけません（エネルギー保存則）
これを数式で表すと以下のようになります。

$$\begin{align*} 1 &\geq\int_{\bold{l}\in\Omega}K(\bold{n}\cdot\bold{l})^+ d\bold{l} \\ 1&\geq \pi K \\ \frac{1}{\pi}&\geq K \end{align*}$$

以上から、エネルギー保存則を満たす完全拡散面の反射率は最大でも$\frac{1}{\pi}$となることがわかります。
材質固有の反射率を$\bold{k}_d$とすると、ランバート反射の反射率$f_d$は

$$f_d(\bold{k}_d)=\frac{1}{\pi}\bold{k}_d$$

以上が正規化ランバートモデルの導出になります。

SmithモデルとGGXモデル

（ちょっと厳密でない表現になりますが）幾何関数Smithモデルは法線分布関数のモデルから算出されます。つまり、法線分布関数にGGXモデルを採用した場合、SmithモデルもGGXに合わせたモデルにする必要があります。実装する際は参考資料をもとに行う場合が多いと思われますが、幾何関数のモデルと法線分布関数のモデルがちゃんと対応しているか確認する必要があります。

このことについて、RealTimeRendering(第4版)を参考に少しだけ詳しく書いてみます。

まず、Smithモデルの$G()$は以下のように表されます

$$G_1(\bold{m},\bold{v})=\frac{\chi^+(\bold{m}\cdot\bold{v})}{1+\Lambda(\bold{v})}$$

$\bold{m}$はマイクロジオメトリの法線ベクトルです。

$\chi^+$は次のような特性関数です。hlslで言うところのstep(0, x)と同じようなものといったらわかりやすいかもしれません。（x==0のときだけ値が違うのに注意してください）

$$\chi^+(x)=\begin{cases} 1, \quad where \; x > 0, \\ 0, \quad where \; x \leq 0, \end{cases}$$

$\Lambda$は法線分布関数によって異なります。GGXの場合は以下。

$$a=\frac{\bold{n}\cdot\bold{s}}{\alpha\sqrt{1-(\bold{n}\cdot\bold{s})^2}}$$

$$\Lambda_{GGX}(a)=\frac{-1+\sqrt{1+\frac{1}{a^2}}}{2}$$

$\bold{s}$は$\bold{l}$か$\bold{v}$のどちらかを表すベクトルです。

参考文献

• Direct3D 12 ゲームグラフィックス実践ガイド
• Physically Based Rendering in Filament
• Background:Physicsand Math of Shading
• リアルタイムレンダリング 第4版
• 【大学数学】立体角(3次元における角度)【解析学】