BRDF理論及shader實現（下）

接上篇：
BRDF理論及shader實現（上）

Specular BRDF

對於specular分量來說，$f_m$是一個遵循菲涅爾反射定律的鏡面BRDF項，此時的$f_m$滿足（[3]和[21]有詳細的推導）：

$$f_m(\mathbf{l},\mathbf{v},\mathbf{m})=F(\mathbf{v},\mathbf{m})\frac{{\delta }_{{\omega }_m}(\mathbf{h},\mathbf{m})}{4(\mathbf{l}\cdot \mathbf{h}{)}^2}$$

$\mathbf{h}$表示half vector，是$\mathbf{v}$和$\mathbf{l}$的平均；這裡分母上第一次出現了$4$，這也是後面specular BRDF公式的分母上的$4$的來源。${\delta }_{{\omega }_m}(\mathbf{s},\mathbf{o})$記為：

$${\int }_{\Omega }g(\mathbf{o}){\delta }_{{\omega }_o}(\mathbf{s},\mathbf{o})d{\omega }_o=\{\begin{array}{ll}g(\mathbf{s}),& \text{if s∈Ω}\\ 0,& \text{otherwise}\end{array}$$

此時，$f_r$可以化簡為為：

$$f_r(\mathbf{l},\mathbf{v})=\frac{D(\mathbf{h},\alpha )G(\mathbf{v},\mathbf{l},\alpha )F(\mathbf{v},\mathbf{h},F_0)}{4(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l})}$$

其中，$\alpha$取決於我們前面提到的粗糙度roughness，具體為

$$\alpha =roughnes{s}^2$$

可以理解為$\alpha$是對粗糙度roughness的一個對映，$\alpha$將多次被用到。

菲涅爾項及$F_0$會在後面詳細介紹，這裡暫時略過。

shader實現：

// #define saturate(x) clamp(x, 0, 1)
// N     = normal;
// V     = normalize(camPos - WorldPos);
// L     = normalize(LightPos - WorldPos));
// H     = normalize(V + L);
// NdotV = saturate(dot(N, V));
// NdotL = saturate(dot(N, L));
// NdotH = saturate(dot(N, H));
// LdotH = saturate(dot(L, H));
// VdotH = saturate(dot(V, H));
vec3 SpecularBRDF(float NdotV, float NdotL, float NdotH, float LdotH, float VdotH, vec3 F0, float roughness) {
    float r   = roughness * roughness;

    float D = Distribution(NdotH, r);
    float V = Geometry(NdotV, NdotL, r);
    vec3  F = Fresnel(VdotH, F0);

    return D * V * F / (4 * NdotV * NdotL);
}

注意，這裡的NdotL，NdotV等都是clamp到0到1的。

接下來具體看$f_r$中每個分量的可能形式都有哪些。

法向分佈函式 D

這一部分介紹3個法向分佈函式的公式，以及一個推廣。

Beckmann

來源[5]，$D_{Beckmann}$假設微表面的法向分佈是以$\mathbf{n}$為均值的高斯分佈，也即$\mathbf{h}$與$\mathbf{n}$越接近，反射的光線越多，$D_{Beckmann}$越大。再結合$D_{Beckmann}$的積分約束，求得：

$$D_{Beckmann}(\mathbf{h},\alpha )=\frac{{\chi }^{+}(\mathbf{n},\mathbf{h})}{\pi {\alpha }^2(\mathbf{n}\cdot \mathbf{h}{)}^4}{e}^{\frac{(\mathbf{n}\cdot \mathbf{h}{)}^2-1}{{\alpha }^2(\mathbf{n}\cdot \mathbf{h}{)}^2}}$$

這個公式有一個很致命的缺陷，那就是當roughness接近於1的時候，$D_{Beckmann}$在$\mathbf{n}\cdot \mathbf{h}\in [0,1]$區間內，不是單調遞增的。表現在渲染上，就是在高光最強的中心點會產生一個暗斑。

shader實現：

float DistributionBeckmann(float NdotH, float r) {
    float NdotH2   = NdotH * NdotH;
    float r2       = r * r;
    float r2NdotH2 = r2 * NdotH2;
    return exp((NdotH2 - 1) / (r2NdotH2)) / (PI * r2NdotH2 * NdotH2);
}

BlinnPhong

來源[6]，BlinnPhong公式純粹是基於經驗的，在恰當選取引數的情況下，它的函式曲線非常接近於Beckmann。

BlinnPhong原始的模型是：

$$D_{Blinn}(\mathbf{h},\alpha )={\chi }^{+}(\mathbf{n},\mathbf{h})\frac{{\alpha }_p+2}{2\pi }(\mathbf{n}\cdot \mathbf{h}{)}^{{\alpha }_p}$$

其中，${\alpha }_p$表示粗糙係數，或者準確的說，是光滑係數——${\alpha }_p$越大，表示物體表面越光滑。

• 當${\alpha }_p=\infty$的時候，表示絕對光滑的物體，此時$D_{Blinn}(\mathbf{h},\alpha )$只有在$\mathbf{h}=\mathbf{n}$，即入射角等於出射角的時候為$\infty$，否則為0。
• 當${\alpha }_p=0$的時候，表示絕對粗糙的物體，$D_{Blinn}(\mathbf{h},\alpha )=\frac{1}{\pi }$，這個式子也是後面會提到的diffuse的式子。

令${\alpha }_p=(\frac{2}{{\alpha }^2}-2)$，則有：

$$D_{Blinn}(\mathbf{h},\alpha )=\frac{{\chi }^{+}(\mathbf{n},\mathbf{h})}{\pi {\alpha }^2}(\mathbf{n}\cdot \mathbf{h}{)}^{(\frac{2}{{\alpha }^2}-2)}$$

這個公式即接近於Beckmann的法向分佈公式，也是常用的BlinnPhong形式。

shader實現：

float DistributionBlinnPhong(float NdotH, float r) {
    float a = r * r;
    return pow(NdotH, 2.0 / a - 2.0) / (PI * a);
}

GGX

來源[3]，GGX是根據實測資料擬合出來的一個公式：

$$D_{GGX}(\mathbf{h},\alpha )=\frac{{\chi }^{+}(\mathbf{n},\mathbf{h})\cdot {\alpha }^2}{\pi ((\mathbf{n}\cdot \mathbf{h}{)}^2({\alpha }^2-1)+1{)}^2}$$

shader實現：

float DistributionGGX(float NdotH, float r) {
    float a2     = r * r;
    float NdotH2 = NdotH * NdotH;
    float nom    = a2;
    float denom  = (NdotH2 * (a2 - 1.0) + 1.0);
    denom        = PI * denom * denom;
    return nom / max(denom, 0.001);
}

除了這三種公式，還有更多更復雜的法向分佈函式D，具體可以參考[17]。但是其實最常用的還是GGX（及其各向異性模式），無論是遊戲還是影視行業都比較喜歡用GGX。

GTR

Burley通過對Berry（與GGX公式類似，分母上的指數為1）和GGX公式的觀察，提出了廣義的Trowbridge-Reitz（Generalized-Trowbridge-Reitz，GTR）法線分佈函式：

$$D_{GTR}(\mathbf{h},\alpha )=\frac{c\cdot {\chi }^{+}(\mathbf{n},\mathbf{h})}{((\mathbf{n}\cdot \mathbf{h}{)}^2({\alpha }^2-1)+1{)}^{\gamma }}$$

其中，$c$表示縮放係數，是一個常數；$\gamma$用於控制尾部的形狀，當$\gamma =1$的時候，$D_{GTR}$就是Berry公式，當$\gamma =2$的時候，$D_{GTR}$就是$D_{GGX}$。

$\gamma$的取值對$D_{GTR}$的影響如下圖所示。

以下是$\gamma =1$和$\gamma =2$時的shader實現：

float DistributionGTR1(float NdotH, float r)
{
    if (r >= 1) return 1/PI;
    float a2 = r*r;
    float t = 1 + (a2-1)*NdotH*NdotH;
    return (a2-1) / (PI*log(a2)*t);
}

float DistributionGTR2(float NdotH, float r)
{
    float a2 = r*r;
    float t = 1 + (a2-1)*NdotH*NdotH;
    return a2 / (PI * t * t);
}

效果對比

可以看出，BlinnPhong和Beckmann的差異不大。而GGX有著更平滑的邊緣和更小的峰值。除此之外，GGX運算壓力更小，因為它沒有指數操作。

遮擋項 G

和法向分佈函式$D$一樣，遮擋項$G$也是入射角、出射角和表面粗糙度的函式。

有些文章會把遮擋項G和BRDF的分母$(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})$放在一起組成一項約分掉，這也是一種優化思路，因為G通常包含這兩個cosine因子。這裡我們約定本文的遮擋項$G$是不約分$(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})$的$G$。

Implicit

來源[7]，有些BRDF公式會忽略遮擋項G，將其跟分母上的$(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})$一起忽略掉，這就有了第一個隱式$G$：

$$G_{Implicit}(\mathbf{l},\mathbf{v},\mathbf{h})=(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})$$

它的形態大概是，當且僅當視角和光源都垂直於物體表面的時候，$G_{Implicit}=1$，光源、視角與物體表面法線的夾角越大，$G_{Implicit}$越小，直到衰減為0，這也是很符合常識的。

shader實現：

float GeometryImplicit(float NdotV, float NdotL) {
    return NdotL * NdotV;
}

但是隱式遮擋項$G_{Implicit}$最大的問題在於，它隨著視角的衰減速度太快，這會使得高光區域太窄。為了解決這個問題，我們繼續看顯式的遮擋項$G$。

Cook-Torrance

來源[9]，$G_{Cook-Torrance}$解決了$G_{Implicit}$衰減速度太快的問題：

$$G_{Cook-Torrance}(\mathbf{l},\mathbf{v},\mathbf{h})=\min \left(1,\frac{2(\mathbf{n}\cdot \mathbf{h})(\mathbf{n}\cdot \mathbf{v})}{\mathbf{v}\cdot \mathbf{h}},\frac{2(\mathbf{n}\cdot \mathbf{h})(\mathbf{n}\cdot \mathbf{l})}{\mathbf{v}\cdot \mathbf{h}}\right)$$

shader實現：

float GeometryCookTorrance(float NdotV, float NdotL, float VdotH, float NdotH) {
    float ct1 = 2 * NdotH * NdotV / VdotH;
    float ct2 = 2 * NdotH * NdotL / VdotH;
    return min(1, min(ct1, ct2));
}

Kelemen

來源[10]，也是解決$G_{Implicit}$衰減速度太快的問題，同時$G_{Kelemen}$比$G_{Cook-Torrance}$的效率更高：

$$G_{Kelemen}(\mathbf{l},\mathbf{v},\mathbf{h})=\frac{(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})}{(\mathbf{v}\cdot \mathbf{h}{)}^2}$$

shader實現：

float GeometryKelemen(float NdotV, float NdotL, float VdotH) {
    return NdotV * NdotL / (VdotH * VdotH);
}

Neumann

來源[8]，$G_{Neumann}$用另一種方式解決了$G_{Implicit}$衰減速度太快的問題：

$$G_{Neumann}(\mathbf{l},\mathbf{v},\mathbf{h})=\frac{(\mathbf{n}\cdot \mathbf{l})(\mathbf{n}\cdot \mathbf{v})}{\max (\mathbf{n}\cdot \mathbf{l},\mathbf{n}\cdot \mathbf{v})}$$

shader實現：

float GeometryNeumann(float NdotV, float NdotL) {
    return (NdotL * NdotV) / max(NdotL, NdotV);
}

但是，以上三個解決方案也不夠完美。前面提到過，遮擋項G應該是入射角、出射角和表面粗糙度的函式，而以上四個G，包括隱式遮擋項都與粗糙度無關。

Smith

Smith家族[13]都是採用了前面介紹的$G_1$相乘的形式：

$$G_2(\mathbf{l},\mathbf{v},\mathbf{h})=G_1(\mathbf{l})G_1(\mathbf{v})$$

他們之間的區別就是$G_1$的選取不同。

Beckmann

Beckmann的$G$是跟$D$一起提出的，前面介紹過$G$是可以從$D$推匯出來的，因此Beckmann的$\Lambda$為：

$$\begin{array}{rl}c& =\frac{\mathbf{n}\cdot \mathbf{v}}{\alpha \sqrt{1-(\mathbf{n}\cdot \mathbf{v}{)}^2}}\\ \Lambda (\mathbf{v})& =\frac{\text{erf}(c)-1}{2}+\frac{1}{2c\sqrt{\pi }}\exp (-c^2)\end{array}$$

但是由於有$\text{erf}$函式的存在，計算起來過於複雜，因此通常用如下的近似形式：

$$\Lambda (\mathbf{v})\approx \{\begin{array}{ll}\frac{1-1.259x+0.396c^2}{3.535c+2.181c^2},& \text{if }c<1.6\\ 0,& \text{if }c\ge 1.6\end{array}$$

因此，Beckmann的$G_1$為

$$G_{Beckmann}(\mathbf{v})\approx \{\begin{array}{ll}\frac{3.535c+2.181c^2}{1+2.276c+2.577c^2},& \text{if }c<1.6\\ 1,& \text{if }c\ge 1.6\end{array}$$

shader實現：

float GeometryBeckmann(float NdotV, float r) {
    float c  = NdotV / (r * sqrt(1 - NdotV * NdotV));
    float c2 = c * c;
    if (c < 1.6)
        return (3.535 * c + 2.181 * c2) / (1 + 2.276 * c + 2.577 * c2);
    else
        return 1.0;
}
float GeometrySmithBeckmann(float NdotV, float NdotL, float r) {
    float ggx2 = GeometryBeckmann(NdotV, r);
    float ggx1 = GeometryBeckmann(NdotL, r);
    return ggx1 * ggx2;
}

GGX

GGX[3]跟Beckmann類似，都是從法向分佈函式推匯出來的：

$$\begin{array}{rl}c& =\frac{\mathbf{n}\cdot \mathbf{v}}{\alpha \sqrt{1-(\mathbf{n}\cdot \mathbf{v}{)}^2}}\\ \Lambda (\mathbf{v})& =\frac{-1+\sqrt{1+\frac{1}{c^2}}}{2}\end{array}$$

對應的$G_1$定義為

$$G_{GGX}(\mathbf{v})=\frac{2(\mathbf{n}\cdot \mathbf{v})}{(\mathbf{n}\cdot \mathbf{v})+\sqrt{{\alpha }^2+(1-{\alpha }^2)(\mathbf{n}\cdot \mathbf{v}{)}^2}}$$

shader實現：

float GeometryGGX(float NdotV, float r) {
    float r2 = r * r;
    return (2 * NdotV) / (NdotV + sqrt(r2 + (1 - r2) * NdotV * NdotV));
}
float GeometrySmithGGX(float NdotV, float NdotL, float r) {
    float ggx2 = GeometryGGX(NdotV, r);
    float ggx1 = GeometryGGX(NdotL, r);
    return ggx1 * ggx2;
}

GGX Joint

前面提到的GGX用的是$G_2=G_1\ast G_1$的separable G，如果用height-correlated G，那麼$G_2$變為：

$$\begin{array}{rl}{G}_{2-GGXJoint}(\mathbf{l},\mathbf{v},\mathbf{m})& =\frac{1}{1+\Lambda (\mathbf{l})+\Lambda (\mathbf{v})}\\ & =\frac{2(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l})}{(\mathbf{n}\cdot \mathbf{l})\cdot \sqrt{{\alpha }^2+(1-{\alpha }^2)(\mathbf{n}\cdot \mathbf{v}{)}^2}+(\mathbf{n}\cdot \mathbf{v})\cdot \sqrt{{\alpha }^2+(1-{\alpha }^2)(\mathbf{n}\cdot \mathbf{l}{)}^2}}\end{array}$$

shader實現：

float GeometrySmithGGXJoint(float NdotV, float NdotL, float r) {
    float r2 = r * r;
    float Vis_SmithV = NdotL * sqrt(NdotV * (NdotV - NdotV * r2) + r2);
	float Vis_SmithL = NdotV * sqrt(NdotL * (NdotL - NdotL * r2) + r2);
	return 2 * NdotV * NdotL / (Vis_SmithV + Vis_SmithL);
}

為了提高計算效率，UE4對GGX Joint方法做了一個近似，公式為：

$$\begin{array}{rl}{G}_{2-GGXJoint}(\mathbf{l},\mathbf{v},\mathbf{m})& =\frac{1}{1+\Lambda (\mathbf{l})+\Lambda (\mathbf{v})}\\ & \approx \frac{2(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l})}{(\mathbf{n}\cdot \mathbf{l})\cdot (\alpha +(1-\alpha )(\mathbf{n}\cdot \mathbf{v}))+(\mathbf{n}\cdot \mathbf{v})\cdot (\alpha +(1-\alpha )(\mathbf{n}\cdot \mathbf{l}))}\end{array}$$

shader實現：

float GeometryGGXJointApprox(float NdotV, float NdotL, float r) {
    return (NdotV) / (NdotL * (r + (1 - r) * NdotV));
}
float GeometrySmithGGXJointApprox(float NdotV, float NdotL, float r) {
	float Vis_SmithV = NdotL * ( NdotV * ( 1 - r ) + r );
	float Vis_SmithL = NdotV * ( NdotL * ( 1 - r ) + r );
	return 2 * NdotV * NdotL / ( Vis_SmithV + Vis_SmithL );
}

Schlick-Beckmann

Schlick[11]的$G_1$定義為
$$\begin{gathered} k=\alpha \sqrt{\frac{2}{\pi }} \\ {G}_{Schlick}(\mathbf{v})=\frac{\mathbf{n}\cdot \mathbf{v}}{(\mathbf{n}\cdot \mathbf{v})(1-k)+k} \end{gathered}$$
shader實現：

float GeometrySchlickBeckmann(float NdotV, float r) {
    float k     = (r)*sqrt(2.0 / PI);
    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;
    return nom / denom;
}
float GeometrySmithSchlickBeckmann(float NdotV, float NdotL, float r) {
    float ggx2 = GeometrySchlickBeckmann(NdotV, r);
    float ggx1 = GeometrySchlickBeckmann(NdotL, r);
    return ggx1 * ggx2;
}

Schlick-GGX

Schlick-GGX[12]曾經是UE4所採用的的一個模型，跟Schlick有些類似，$G_1$定義為
$$\begin{gathered} k=\frac{\alpha }{2} \\ {G}_{Schlick}(\mathbf{v})=\frac{\mathbf{n}\cdot \mathbf{v}}{(\mathbf{n}\cdot \mathbf{v})(1-k)+k} \end{gathered}$$
shader實現：

float GeometrySchlickGGX(float NdotV, float r) {
    float k     = r * 0.5;
    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;
    return nom / denom;
}
float GeometrySmithSchlickGGX(float NdotV, float NdotL, float r) {
    float ggx2 = GeometrySchlickGGX(NdotV, r);
    float ggx1 = GeometrySchlickGGX(NdotL, r);
    return ggx1 * ggx2;
}

這裡面還有一個細節，那就是迪士尼後來提出了對粗糙粗roughness做一個remapping，使得它更接近於真實：
$${\alpha }^{\prime }=(\frac{roughness+1}{2}{)}^2$$
其他的部分不變。這樣shader實現為：

float GeometrySmithSchlickGGX(float NdotV, float NdotL, float roughness) {
    float r = (roughness + 1.0) * 0.5; // remapping roughness
    r = r * r
    float ggx2 = GeometrySchlickGGX(NdotV, r);
    float ggx1 = GeometrySchlickGGX(NdotL, r);
    return ggx1 * ggx2;
}

注意，此時GeometrySmithSchlickGGX的輸入引數不是r，而改為了roughness。

優化

考慮到幾乎所有$G$都帶有$(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l})$項，可以跟$f_r(\mathbf{l},\mathbf{v})$的分母約分，因此在實現時，可以考慮定義
$$G^{\prime }=\frac{G}{(\mathbf{n}\cdot \mathbf{v})(\mathbf{n}\cdot \mathbf{l})}$$
節省一部分計算。

這樣做不只是出於效能的考慮，也是出於精度的考慮。如果$\mathbf{n}\cdot \mathbf{v}$和$\mathbf{n}\cdot \mathbf{l}$的乘積接近於0，那麼specular項的分母會非常小，嚴重影響其精度，極端的情況下會在過渡區域產生一道割裂的分界線。下圖展示了$(\mathbf{n}\cdot \mathbf{v})\ast (\mathbf{n}\cdot \mathbf{l})\ast 10$（左）、未優化時（中）、優化後（右）的效果，可以看出左側兩張圖的分界線非常吻合。優化後則沒有顏色割裂的問題。

效果對比

取roughness = 0.9，計算球體的遮擋項效果為：

最後一排是常用的幾種方法，差異並不大，邊緣的過渡也比較好。

菲涅爾項 F

菲涅爾項描述的是物體表面的反射、折射現象。一般我們會採用常量$F_0$來計算菲涅爾項$F(\mathbf{v},\mathbf{h},F_0)$。

要說明白菲涅爾項，得從光學在介質表面的折射反射現象說起。我們知道光線會在介質表面產生不連續性，具體表現為一部分光線反射——遵循光線反射定律，入射角等於反射角；另一部分光線會折射進入介質——遵循光線折射定律，折射角取決於入射角的大小以及構成交介面的兩種材質，即斯涅耳定律（Snell’s law）：

$$n_1\sin ({\theta }_i)=n_2\sin ({\theta }_t)$$

斯涅耳定律描述的僅僅是光線的角度，但是圖形學研究的其實是光線的radiance/irradiance，所以我們要更進一步。定義Fresnel reflectance $R_F$為反射光線的radiance佔入射光線radiance的比例，$R_F$是入射角${\theta }_i$的函式。那麼對於入射光線$L_i$，在角度${\theta }_i$時反射光線的radiance為$R_F({\theta }_i)L_i$。再考慮折射部分，根據能量守恆，沒有反射的能量都會被折射（不考慮被吸收的能量），因此折射的flux佔入射flux的比例是$1-R_F$。這裡需要強調的是，radiance定義的是“irradiance每立體角”，它的大小跟角度有關係，因此折射光線的radiance $L_t$不能簡單用$1-R_F$乘上$L_i$，而要轉換角度：
$$L_t=(1-R_F({\theta }_i))\frac{{\sin }^2{\theta }_i}{{\sin }^2{\theta }_t}{L}_i$$
將斯涅耳定律帶入上式，得到：
$$L_t=(1-R_F({\theta }_i))\frac{n_2^2}{n_1^2}{L}_i$$
介紹了這麼多$R_F$的相關知識，其實關鍵點還是前面說的，$R_F$是入射角${\theta }_i$的函式。我們再回頭考慮這個$R_F$與入射角${\theta }_i$的關係。當${\theta }_i=90°$的時候，即$R_F(90°)$，此時入射光平行於平面，垂直於法向，不存在折射光線，$R_F(90°)=1$；當${\theta }_i=0°$的時候，即$R_F(0°)$，此時反射光線佔比最低，根據不同的材質這個$R_F(0°)$有不同的值，Real-time Rendering[14]給出了常見的材質的$R_F$與${\theta }_i$的關係曲線：

為了近似這個曲線，採取的策略是利用$R_F(0°)$，也就是前面說的$F_0$：
$$R_F({\theta }_i)\approx R_F(0°)+(1-R_F(0°))(1-\cos {\theta }_i{)}^5$$
這裡有一個預設的假設是，$R_F(90°)=1$，如果，$R_F(90°)$未知，$R_F({\theta }_i)$應該寫為：
$$R_F({\theta }_i)\approx R_F(0°)+(R_F(90°)-R_F(0°))(1-\cos {\theta }_i{)}^5$$
這個$R_F(90°)$也就是$F_{90}$。

最後，我們看一下$F_0$怎麼計算。對於dielectrics來說，$F_0$的值取決於折射率，公式為：
$$F_0=0.16\cdot reflectanc{e}^2$$
其中，$reflectance$由物體表面的材質定義。

對於dielectric，$F_0$通過金屬度metallic和basecolor來計算：
$$F_0=baseColor\cdot metallic$$
綜合dielectrics和dielectric，得到：

    vec3 F0 = 0.16 * reflectance * reflectance * (1.0 - metallic) + baseColor.xyz * metallic;

說明白了$F_0$，我們接下來看看菲涅爾函式$F$有哪些形式。

簡單形式

最簡單的情況，直接令菲涅爾函式等於$F_0$：
$$F_{None}(\mathbf{v},\mathbf{h})=F_0$$
shader實現：

vec3 Fresnel(vec3 F0) {
    return F0;
}

Schlick

來源[11]，公式：
$$F_{Schlick}(\mathbf{v},\mathbf{h})=F_0+(1-F_0)(1-(\mathbf{v}\cdot \mathbf{h}){)}^5$$
也就是我們前面說到的對$R_F$的擬合。shader實現：

vec3 FresnelSchlick(float VdotH, vec3 F0) {
    return F0 + (1.0 - F0) * pow(1.0 - VdotH, 5.0);
}

如果引入$F_{90}$，則變成：
$$F_{Schlick}(\mathbf{v},\mathbf{h})=F_0+(F_{90}-F_0)(1-(\mathbf{v}\cdot \mathbf{h}){)}^5$$
shader實現：

vec3 FresnelSchlick(float VdotH, vec3 F0, vec F90) {
    return F0 + (F90 - F0) * pow(1.0 - VdotH, 5.0);
}

對specular來說，$F_{90}$可以從$F_0$計算得來[1]：

    float F90 = saturate(dot(F0, vec3(50.0 * 0.33)));

Cook-Torrance

來源[9]，公式：
$$\begin{array}{rl}\eta & =\frac{1+\sqrt{F_0}}{1-\sqrt{F_0}}\\ c& =\mathbf{v}\cdot \mathbf{h}\\ g& =\sqrt{{\eta }^2+c^2-1}\\ F_{Cook-Torrance}(\mathbf{v},\mathbf{h})& =\frac{1}{2}{\left(\frac{g-c}{g+c}\right)}^2\left(1+{\left(\frac{(g+c)c-1}{(g-c)c+1}\right)}^2\right)\end{array}$$
shader實現：

float FresnelCookTorrance(float VdotH, float F0) {
    float sqrtF = sqrt(F0);
    float Eta   = (1.0 + sqrtF) / (1.0 - sqrtF);
    float g     = sqrt(Eta * Eta + VdotH * VdotH - 1.0);
    return 0.5 * pow((g - VdotH) / (g + VdotH), 2) *
           (1 + pow(((g + VdotH) * VdotH - 1.0) / ((g - VdotH) * VdotH + 1.0), 2));
}

Diffuse BRDF

相比於繁瑣的specular部分，diffuse部分就簡單的多。diffuse部分由baseColor和diffuse係數相乘得到，即：
$$L_d(\mathbf{v})={\mathbf{c}}_{diff}\cdot f_d$$
shader實現：

    vec3 colorDiffuse = baseColor * DiffuseBRDF(NdotV, NdotL, LdotH, roughness);

接下來看一下$f_d$的可能取值。

Lambert

Lambert模型認為既然diffuse是漫反射，不如簡單地認為各個方向都是一樣的值，即出射光線的radiance與入射光線的角度無關。

$$f_d=\frac{1}{\pi }$$

這個實現相當於，採用BlinnPhong的法向分佈$D_{Blinn}(\mathbf{h},\alpha )=\frac{1}{\pi }$，同時令遮擋項為隱式形式，並且菲涅爾項為1。雖然簡單，但是已經足夠近似現實了，效果還不錯。

shader實現：

float DiffuseLambert() {
    return 1.0 / PI;
}

雖然Lambert模型已經足夠接近真實情況，但是它還是不夠理想。我們前面提到過，diffuse分量本質上是光線折射進入物體表面，經過多次反射再折射出來的現象，也就是它不是物理上真實存在的一個光學現象。而在討論specular菲涅爾項的時候又提到過，反射部分會隨著入射光線的角度變化，那麼折射部分相應的也會隨著入射角度變化，既然如此，來自於折射部分的diffuse分量肯定也是會隨著入射光線的角度而改變的！也就是說，$f_d$是入射角$\mathbf{l}$的函式：$f_d(\mathbf{l})$。

同時，$f_d$也應該是出射角$\mathbf{v}$的函式[14]：$f_d(\mathbf{l},\mathbf{v})$。因為菲涅爾項考慮的是鏡面反射，入射角等於出射角，而diffuse項的入射角不一定等於出射角，因此兩個角度都會影響$f_d$。

再者，前面影響specular分量的引數當中，$roughness$也會影響$f_d$。根據常識，不同粗糙程度的物體的diffuse是有明顯的不同的。即$f_d(\mathbf{l},\mathbf{v},roughness)$。

基於這一點洞察，又有一些新的diffuse模型被提出，希望解決Lambert模型的不足。。

Oren–Nayar

Oren-Nayar模型是對Lambert模型的推廣。[18]指出，Lambert模型對於光滑物體或許還成立，但是對於粗糙物體是不正確的。粗糙的物體在光照下會顯得很平坦，而Lambert模型沒有表現出這種平坦。為了達到這個效果，Oren-Nayar加強了掠射逆反射（入射角和出射角在幾乎同一個方向，並且垂直於法向的情形）的強度。

Oren-Nayar公式如下。

$$\begin{array}{rl}{f}_d& =\frac{1}{\pi }\cdot (A+B\cdot \max (0,\cos \varphi )\cdot \sin \alpha \cdot \tan \beta )\\ A& =1.0-0.5\frac{\alpha }{\alpha +0.33}\\ B& =0.45\frac{\alpha }{\alpha +0.09}\\ \alpha & =\max (\mathbf{l}\cdot \mathbf{n},\mathbf{v}\cdot \mathbf{n})\\ \beta & =\min (\mathbf{l}\cdot \mathbf{n},\mathbf{v}\cdot \mathbf{n})\end{array}$$

其中，$\varphi$表示$\mathbf{l}\mathbf{n}$平面和$\mathbf{v}\mathbf{n}$的夾角。

可以看出，當$roughness-0$的時候，$A=1,B=0$，此時Oren-Nayar模型退化為Lambert模型。

下圖[18]展示了真實照片、Lambert模型與Oren-Nayar的對比。

Hanrahan-Krueger

Hanrahan-Krueger模型[19]其實是源自次表面散射理論，是用於表現次表面散射現象的一個模型。它跟Oren-Nayar模型一樣，對掠射角進行了補償。但是它的補償過於平坦，沒有給出足夠強的峰值，也不太完美。

Hanrahan-Krueger模型和Oren-Nayar模型都不太常用，因此不再贅述。

Burley

Oren–Nayar模型雖然提高了粗糙物體的真實性，但是它對掠射逆反射現象的修正還是不夠真實。為了研究真實材料的物理特性，我們需要一個材質資料庫。

MERL BRDF Database就是這樣一個資料庫。它是由MERL(Mitsubishi Electric Research Laboratories)實驗室建立了的，測量並記錄了不同角度的光源、觀測視角情況下的BRDF數值，考慮到各向異性，每個材質都取樣了90（光源）*90（視角）*180（各向異性）三個維度的資料。如果只考慮各向同性材質，可以將BRDF資料壓縮到一張圖片裡。

如上圖所示，橫軸${\theta }_h$表示half vector $\mathbf{h}$與法向量$\mathbf{n}$之間的夾角。縱軸表示入射角與$\mathbf{h}$的夾角。

Disney通過分析MERL BRDF Database，提出了兩個Lambert模型與事實不符的地方：

1. diffuse也會有類似於specular的光斑；
2. 部分材質的diffuse會在掠射角有明顯的光環，這個現象即掠射逆反射（grazing retroreflection）；

為了解決這些問題，Disney提出了一個diffuse BRDF公式[15]：

$$\begin{array}{rl}{f}_d(\mathbf{l},\mathbf{v})& =\frac{1}{\pi }{F}_{Schlick}(\mathbf{n},\mathbf{l},1,f_{90})F_{Schlick}(\mathbf{n},\mathbf{v},1,f_{90})\\ F_{Schlick}(\mathbf{n},\mathbf{l},f_0,f_{90})& =F_0+(F_{90}-F_0)(1-(\mathbf{n}\cdot \mathbf{l}){)}^5\\ f_{90}& =0.5+2\cdot roughness\cdot {\cos }^2({\theta }_d)\end{array}$$

其中${\theta }_d$是光線$L$和half vector $h$的夾角。這個公式考慮到了入射角和出射角以及粗糙度，並且用類似菲涅爾項的公式（cosine項的五次方）來擬合衰減情況。

shader實現:

float FresnelSchlick(float VdotH, float F0, float F90) {
    return F0 + (F90 - F0) * pow(1.0 - VdotH, 5.0);
}
float DiffuseBurley(float NdotV, float NdotL, float LdotH, float roughness) {
    float f90          = 0.5 + 2.0 * roughness * LdotH * LdotH;
    float lightScatter = FresnelSchlick(NdotL, 1.0, f90);
    float viewScatter  = FresnelSchlick(NdotV, 1.0, f90);
    return lightScatter * viewScatter * (1.0 / PI);
}

總結

BRDF作為渲染裡邊最基礎的知識點，發展的已經相對成熟，雖然偶爾也會有一些改進，但是基本上都是在效率與效能之間做權衡。對於基本BRDF公式的選擇，UE4和Disney有著各自不同的邏輯：

$$\begin{array}{ccccc}& \text{Diffuse BRDF}& \text{Distribution}& \text{Visibility}& \text{Fresnel}\\ \text{UE4}& Lambert& \text{GGX}& GGXJoint(Approx)& \text{Schlick}\\ \text{Disney}& Burley& \text{GGX}& GGX& \text{Schlick}\end{array}$$

斜體表示二者不同的部分。可以看出，UE4選擇的都是高效的模型，而Disney選擇的都是複雜而準確的模型。

個人理解這些差異都是源於UE4和disney應用場景的不同，UE4希望每個模型儘可能高效，因此會拆分開來，針對性優化，比如它單獨設計了針對眼睛的Eye模型，專門渲染毛髮的Hair模型，專門渲染皮膚的subsurface模型等等。而Disney的訴求在於模型的表達力要足夠強，效率反而不那麼重要。

未涉及話題…

本文主要集中在BRDF項的各種實現，順帶介紹了BRDF和微表面理論。還有一些與之相關或更深入，但是沒有涉及到的方向，例如

1. 輻射度量學基礎；
2. BSDF，BTDF等BRDF的進階模型；
3. 各向異性BRDF，subsurface、clearCoat等模型；
4. 環境光、全域性光照等；

篇幅問題，這些方向也無法展開。行文至此，強推圖形學屆的武林祕籍的目錄——Real-time Rendering，此書目前已經出到第四版了，文末也給出了電子書連結[14]。本文涉及的話題書中都有比較深入、全面的介紹。即使RTR不能滿足你，它還提供了多達1000+篇的參考文獻供學習，畢竟“目錄”，名副其實。

參考資料

1. Filament文件，Filament是一個Google寫的用在Android上的PBR渲染器，它的文件非常完善，特別每個BRDF的理論和實現。同時也可以參考它的原始碼，對照學習。
2. Specular BRDF Reference：這個部落格列出了幾大主流specular BRDF的公式，可以作為參考。
3. Walter et al. 2007, Microfacet models for refraction through rough surfaces
4. LearningOpenGL: PBR Theory：這也是一個不錯的學習PBR的教材，有一個PBR的OpenGL實現，以及簡單的理論介紹。
5. Beckmann 1963, The scattering of electromagnetic waves from rough surfaces
6. Blinn 1977, Models of light reflection for computer synthesized pictures
7. Hoffman 2013, Background: Physics and Math of Shading
8. Neumann et al. 1999, Compact metallic reflectance models
9. Cook and Torrance 1982, A Reflectance Model for Computer Graphics
10. Kelemen 2001, A microfacet based coupled specular-matte brdf model with importance sampling
11. Schlick 1994, An Inexpensive BRDF Model for Physically-Based Rendering
12. Karis 2013, Real Shading in Unreal Engine 4
13. Smith 1967, Geometrical shadowing of a random rough surface
14. Real-time Rendering, 4th edition，需要說明的一點是，此書的第四版比第三版增加了很多對BRDF公式的推導和歷史介紹，更具有參考價值。
15. Brent Burley. 2012. Physically Based Shading at Disney. Physically Based Shading in Film and Game Production, ACM SIGGRAPH 2012 Courses.
16. Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs
17. SIGGRAPH 2013 Course, Background: Physics and Math of Shading
18. Generalization of Lambert’s reflectance model
19. Reflection from Layered Surfaces due to Subsurface Scattering
20. SIGGRAPH 2013 Course, Physically Based Shading at Disney
21. PBR Diffuse Lighting for GGX+Smith Microsurfaces