GLSLで物理ベースシェーディング


GLSLで物理ベースシェーディングを試すために、いろいろなBRDFを実装してみました。
![mix.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f7be0e2a-3457-18e2-bf2c-fd291d182299.png)

サンプルを以下においておきました。マクロで使用するBRDFを変更することができます。
http://glslsandbox.com/e#54592.0

マウス位置のx方向で`roughness`パラメータを、y方向で`metallic`パラメータを変更することができます。

## 下準備

BRDFを定義する際に使用する変数です。

```math
\begin{eqnarray}
\rho_{d}&:& 拡散リフレクタンス\\
\rho_{s}&:& 鏡面リフレクタンス\\
\vec{n}&:& 法線(normal)ベクトル\\
\vec{v}&:& 視線(view)ベクトル\\
\vec{l}&:& 光源(light)ベクトル\\
\vec{r}&:& 反射(reflect)ベクトル\\
\vec{h}&:& ハーフ(half)ベクトル
\end{eqnarray}
```

ハーフベクトルは2つのベクトルの中間のベクトルです。
シェーディングの文脈では視線ベクトル$\vec{v}$と光源ベクトル$\vec{l}$の中間ベクトルになり、以下の計算で求めます。

```math
\vec{h} = \frac{\vec{v}+\vec{l}}{|\vec{v}+\vec{l}|}
```

GLSLでベクトルを定義する箇所は以下のようになります。シェーディングでは2つのベクトルがなす角度が必要となることが多いのであらかじめ内積を計算しています。

```glsl
vec3 normal; // 法線ベクトル
vec3 viewDir; // 視線ベクトル
vec3 lightDir; // 光源ベクトル
vec3 reflectDir = reflect(-viewDir, normal); // 視線ベクトル
vec3 halfDir = normalize(viewDir + lightDir); // ハーフベクトル

float dotNV = saturate(dot(normal, viewDir));
float dotNL = saturate(dot(normal, lightDir));
float dotNH = saturate(dot(normal, halfDir));
float dotLH = saturate(dot(lightDir, halfDir));
float dotRL = saturate(dot(reflectDir, lightDir));
```

## 拡散反射BRDF

サンプルでは、以下のマクロで適用する拡散反射BRDFを選択することができます。

```glsl
// 0: normalized lambert
// 1: difuse disney
// 2: normalized diffuse disney
// 3: oren nayar
// others: no diffuse
#define DIFFUSE_BRDF 0
```

拡散反射BRDFの影響のみを見たい場合は以下のように設定して、鏡面反射成分をオフにしてください。

```glsl
#define SPECULAR_BRDF -1
```

### 正規化Lambert

![lambert.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/fb0ccbb1-d4e9-905a-d36b-bd395e594ca8.png)

Lambert拡散反射をエネルギー保存則を守るように正規化したものです。BRDFは以下のようになります。

```math
f_{d} = \frac{\rho_d}{\pi}
```

```glsl
vec3 diffuseMormalizedLambertBrdf(vec3 reflectance) {
    return reflectance * INV_PI;
}
```

### Disney Diffuse

![disney.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f7311547-585f-17c5-96eb-47af43cad5a2.png)

フレネル反射率を考慮した拡散反射BRDFです。$roughness$は表面の粗さを表すパラメータです。BRDFは視線方向のフレネル反射率とライト方向のフレネル反射率を乗算したものとなっています。

```math
f_{D90} = 0.5+2(\vec{h}\cdot\vec{l})^{2}roughness\\
f_{d} = \frac{\rho_d}{\pi}(1+(f_{D90}-1)(1-\vec{n}\cdot\vec{l})^5)(1+(f_{D90}-1)(1-\vec{n}\cdot\vec{v})^5)
```

GLSLで実装すると以下のようになります。

```glsl
float fresnelSchlick(float f0, float f90, float cosine) {
    return f0 + (f90 - f0) * pow(1.0 - cosine, 5.0);
}

vec3 diffuseDisneyBrdf(vec3 reflectance, float dotNV, float dotNL, float dotLH, float roughness) {
    float fd90 = 0.5 + 2.0 * dotLH * dotLH * roughness;
    float fl = fresnelSchlick(1.0, fd90, dotNL);
    float fv = fresnelSchlick(1.0, fd90, dotNV);
    return reflectance * INV_PI * fl * fv;
}
```

**参考**

* [拡散BRDF \- Qiita](https://qiita.com/mebiusbox2/items/1cd65993ffb546822213)
* [物理ベースレンダリングを柔らかく説明してみる（２） \- Qiita](https://qiita.com/emadurandal/items/ca6406fb20f4878ab222)

### Normalized Disney Diffuse

![normalized-disney.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/144af923-be72-ef30-4ca6-4ef3440b1cda.png)

Disney Diffuse BRDFは$roughness$が大きくグレージング角に近づくと、値が1を超えてエネルギー保存則を守らなくなります。これを防ぐために正規化したものが以下になります。

```glsl
vec3 diffuseNormalizedDisneyBrdf(vec3 reflectance, float dotNV, float dotNL, float dotLH, float roughness) {
    float bias = mix(0.0, 0.5, roughness);
    float factor = mix(1.0, 1.0 / 1.51, roughness);
    float fd90 = bias + 2.0 * dotLH * dotLH * roughness;
    float fl = fresnelSchlick(1.0, fd90, dotNL);
    float fv = fresnelSchlick(1.0, fd90, dotNV);
    return reflectance * INV_PI * fl * fv * factor;
}
```

**参考**

* [拡散BRDF \- Qiita](https://qiita.com/mebiusbox2/items/1cd65993ffb546822213)

### Oren Nayar

![oren-nayar.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/7cdaa892-9f95-af70-7244-3f74f7394d5e.png)

マイクロファセットを考慮した拡散反射BRDFです。$roughness$は表面の粗さを表すパラメータです。

```math
f_{d}=\frac{\rho_d}{\pi}(A+B\cdot max(0, cos(\phi_l-\phi_v))\cdot sin\alpha\cdot tan\beta)\\
A=1-0.5\frac{roughness^2}{roughness^2+0.33}\\
B=0.45\frac{roughness^2}{roughness^2+0.09}\\
\alpha=max(\theta_{l}, \theta{v})\\
\beta=min(\theta_{l}, \theta{v})
```

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

```glsl
vec3 diffuseOrenNayarBrdf(vec3 reflectance, vec3 normal, vec3 viewDir, vec3 lightDir, float roughness) {
    float dotNV = saturate(dot(normal, viewDir));
    float dotNL = saturate(dot(normal, lightDir));
    float roughness2 = roughness * roughness;
    float a = 1.0 - 0.5 * roughness2 / (roughness2 + 0.33);
    float b = 0.45 * roughness2 / (roughness2 + 0.09);
    float cosPhi = dot(normalize(viewDir - dotNV * normal), normalize(lightDir - dotNL * normal)); // cos(phi_v, phi_l)
    float sinNV = sqrt(1.0 - dotNV * dotNV);
    float sinNL = sqrt(1.0 - dotNL * dotNL);
    float s = dotNV < dotNL ? sinNV : sinNL; // sin(max(theta_v, theta_l))
    float t = dotNV > dotNL ? sinNV / dotNV : sinNL / dotNL; // tan(min(theta_v, theta_l))
    return reflectance * INV_PI * (a + b * cosPhi * s * t);
}
```

**参考**

* [オーレン・ネイヤー反射 \- Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%83%BC%E3%83%AC%E3%83%B3%E3%83%BB%E3%83%8D%E3%82%A4%E3%83%A4%E3%83%BC%E5%8F%8D%E5%B0%84)
* [t\-pot『粘土ライティング』](https://t-pot.com/program/80_OrenNayer/index.html)

## 鏡面反射BRDF

サンプルでは、以下のマクロで適用する鏡面反射BRDFを選択することができます。

```glsl
// 0: noramlized phong
// 1: normalized blinn-phong
// 2: cook-torrance
// others: no specular
#define SPECULAR_BRDF 0
```

鏡面反射BRDFの影響のみを見たい場合は以下のように設定して、拡散反射成分をオフにしてください。

```glsl
#define DIFFUSE_BRDF -1
```

### 正規化Phong

![phong.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/b018801d-2601-eaff-1d8a-b3dbb891e21a.png)

Phong鏡面反射をエネルギー保存則を守るように正規化したBRDFです。通常のPhong反射のBRDFに正規化項$(n + 1) / 2\pi$がついた形になっています。

```math
f_{s} = \rho_{s}(\vec{r}\cdot\vec{l})^n\frac{n + 1}{2\pi}
```

GLSLで実装すると以下のようになります。

```glsl
vec3 specularNormalizedPhongBrdf(vec3 reflectance, float dotRL, float power) {
    float norm = (power + 1.0) * INV_TWO_PI;
    return reflectance * pow(dotRL, power) * norm;
}
```

**参考**

* [☆PROJECT ASURA☆ \[Direct3D 11\] 『正規化Phong』](http://www.project-asura.com/program/d3d11/d3d11_006.html)


### 正規化Blinn-Phong

![blinn-phong.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/2621a9a9-57b5-8a8b-d59a-b0b4a23dfe81.png)


Blinn-Phong鏡面反射をエネルギー保存則を守るように正規化したBRDFです。通常のBlinn-Phong反射のBRDFに正規化項$(n + 2)/2\pi$がついた形になっています。

```math
f_{s} = \rho_{s}(\vec{n}\cdot\vec{h})^n\frac{n + 2}{2\pi}
```

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

```glsl
vec3 specularNormalizedBlinnPhongBrdf(vec3 reflectance, float dotNH, float power) {
    float norm = (power + 2.0) * INV_TWO_PI;
    return reflectance * pow(dotNH, power) * norm;
}
```

**参考**

* [Blinn\-Phong スペキュラの正規化について \- OLD hanecci’s blog : 旧 はねっちブログ](http://d.hatena.ne.jp/hanecci/20130505/p2)

### Cook Torrance

![cook-torrance.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/01f9caa7-6bed-9a80-32c3-0f058433189c.png)

マイクロファセットを考慮した鏡面反射BRDFです。$F$がフレネル項、$D$が法線分布関数、$G$が幾何減衰項です。

```math
f_{s}=\frac{FDG}{4(\vec{n}\cdot\vec{v})(\vec{n}\cdot\vec{l})}
```

#### フレネル項

フレネル項はどれぐらい入射した光が鏡面反射するかを表しています。ここではShlickによる近似式を用います。$F0$は垂直方向から入射した光の反射率を表しており、グレージング角に近づくほど反射率が大きくなります。

```math
F(\vec{v}, \vec{h}) = F_{0} + (1 - F_{0})(1 - \vec{v}\cdot\vec{h})^5 
```

```glsl
vec3 fresnelSchlick(vec3 f0, float cosine) {
    return f0 + (1.0 - f0) * pow(1.0 - cosine, 5.0);
}
```

#### 法線分布関数

法線分布関数はマイクロファセットのミクロな面がどのように分布しているかを表しています。
ここではGGXという法線分布関数を利用しています。

```math
D(\vec{h})=\frac{roughness^2}{\pi(1-(1-roughness^2)(\vec{n}\cdot\vec{h})^2)^2}
```

```glsl
float normalDistributionGgx(vec3 normal, vec3 halfDir, float roughness) {
    float roughness2 = roughness * roughness;
    float dotNH = saturate(dot(normal, halfDir));
    float a = (1.0 - (1.0 - roughness2) * dotNH * dotNH);
    return roughness2 * INV_PI / (a * a);
}
```

#### 幾何減衰項項

幾何減衰項は入射する光がマイクロファセットの幾何構造でどれぐらい遮蔽(Masking&Shading)されるかを表しています。
ここでは、Smith Masking Shadowing FunctionとSmith Joint Masking Shadowing Functionの2つを実装しています。

##### Smith Masking Shading Function

```math
G(\vec{v},\vec{l})=\frac{1}{1+\Lambda(\vec{v})}\cdot\frac{1}{1+\Lambda(\vec{l})}\\
\Lambda(\vec{x})=\frac{-1+\sqrt{}1+roughness^2(\frac{1}{(\vec{x}\cdot\vec{n})^2}-1)}{2}
```

```glsl
float maskingShadowingSmith(vec3 normal, vec3 viewDir, vec3 lightDir, float roughness) {
    float roughness2 = roughness * roughness;
    float dotNV = saturate(dot(normal, viewDir));
    float dotNL = saturate(dot(normal, lightDir));
    float lv = 0.5 * (-1.0 + sqrt(1.0 + roughness2 * (1.0 / (dotNV * dotNV) - 1.0)));
    float ll = 0.5 * (-1.0 + sqrt(1.0 + roughness2 * (1.0 / (dotNL * dotNL) - 1.0)));
    return (1.0 / (1.0 + lv)) * (1.0 / (1.0 + ll));
}
```

##### Smith Joint Masling Shading Function

```math
G(\vec{v},\vec{l})=\frac{1}{1+\Lambda(\vec{v})+\Lambda(\vec{l})}
```

```glsl
float maskingShadowingSmithJoint(vec3 normal, vec3 viewDir, vec3 lightDir, float roughness) {
    float roughness2 = roughness * roughness;
    float dotNV = saturate(dot(normal, viewDir));
    float dotNL = saturate(dot(normal, lightDir));
    float lv = 0.5 * (-1.0 + sqrt(1.0 + roughness2 * (1.0 / (dotNV * dotNV) - 1.0)));
    float ll = 0.5 * (-1.0 + sqrt(1.0 + roughness2 * (1.0 / (dotNL * dotNL) - 1.0)));
    return 1.0 / (1.0 + lv + ll);
}
```

---

Cook Torrance鏡面反射BRDFのGLSL実装は以下のようになります。


```glsl
vec3 specularCookTorranceBrdf(vec3 reflectance, vec3 normal, vec3 viewDir, vec3 lightDir, float roughness) {
    vec3 halfDir = normalize(viewDir + lightDir);
    float dotNV = saturate(dot(normal, viewDir));
    float dotNL = saturate(dot(normal, lightDir));
    float dotVH = saturate(dot(viewDir, halfDir));
    float d = normalDistributionGgx(normal, halfDir, roughness);
    #if COOK_TORRANCE_MASKING_SHADOWING_FUNCITON == 0
    float g = maskingShadowingSmith(normal, viewDir, lightDir, roughness);
    #else 
    float g = maskingShadowingSmithJoint(normal, viewDir, lightDir, roughness);
    #endif
    vec3 f = fresnelSchlick(reflectance, dotVH);
    return d  * g * f / (4.0 * dotNV * dotNL);
}
```

以下のマクロで使用する幾何減衰項を変更することできます。

```
// 0: smith masking shadowing function
// 1: smith joint masking shadowing function
#define COOK_TORRANCE_MASKING_SHADOWING_FUNCITON 0
```

使用する法線分布関数によってTorrance Sparrow鏡面反射BRFと呼んだりしているようですが、よくわかりませんでした...

**参考**

* [脱・完全鏡面反射～GGXについて調べてみた～ \- Qiita](https://qiita.com/_Pheema_/items/f1ffb2e38cc766e6e668)

## RoughnessとSmoothnessについて

表面の粗さを表す`roughness`というパラメータを使用するBRDFがありました。一方でライブラリや3Dエンジンではパラメータを`smoothness`としている場合もありますが、基本的な考え方は`roughness=1-smoothness`でよさそうです。またパラメータに対して見た目がリニアに変化するように、実際の計算では`roughness^2`を`roughness`として使用する場合もあるようです。

## 拡散反射BRDF + 鏡面反射BRDF

拡散反射BRDFと鏡面反射BRDFを以下のように組み合わせてBRDFを定義します。$k$には鏡面反射の割合を表すフレネル反射率を使うのが物理的に正しそうですが、`metalllic`というパラメータで任意に制御できるようにしている場合が多いような気がします。

```math
f = (1.0 - k)f_{d} + kf_{s}
```

実際には非金属でも4%は反射するので、Unityではそのあたりも考慮して計算されています。
また、Cook Torrance BRDFのようにBRDFでフレネル係数を考慮している場合があるので、そこもうまくやらないと駄目そうです。よくわかってないですが...

**参考**

* [基礎からはじめる物理ベースレンダリング（実装編） \- Qiita](https://qiita.com/mebiusbox2/items/8a4734ab5b0854528789)

## 参考

* [\[Unity\] PBR（物理ベースレンダリング）をやってみる \- Qiita](https://qiita.com/edo_m18/items/636ee3e54e8cc72b1241)
* [拡散BRDF \- Qiita](https://qiita.com/mebiusbox2/items/1cd65993ffb546822213)
* [物理ベースレンダリングを柔らかく説明してみる（２） \- Qiita](https://qiita.com/emadurandal/items/ca6406fb20f4878ab222)
* [オーレン・ネイヤー反射 \- Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%83%BC%E3%83%AC%E3%83%B3%E3%83%BB%E3%83%8D%E3%82%A4%E3%83%A4%E3%83%BC%E5%8F%8D%E5%B0%84)
* [t\-pot『粘土ライティング』](https://t-pot.com/program/80_OrenNayer/index.html)
* [☆PROJECT ASURA☆ \[Direct3D 11\] 『正規化Phong』](http://www.project-asura.com/program/d3d11/d3d11_006.html)
* [Blinn\-Phong スペキュラの正規化について \- OLD hanecci’s blog : 旧 はねっちブログ](http://d.hatena.ne.jp/hanecci/20130505/p2)
* [脱・完全鏡面反射～GGXについて調べてみた～ \- Qiita](https://qiita.com/_Pheema_/items/f1ffb2e38cc766e6e668)
* [基礎からはじめる物理ベースレンダリング（実装編） \- Qiita](https://qiita.com/mebiusbox2/items/8a4734ab5b0854528789)