LATEX

Rendering Equaltion

$$ \def\red#1{\textcolor[rgb]{1,0,0}{#1}} \def\blue#1{\textcolor[rgb]{0,1,0}{#1}} \def\green#1{\textcolor[rgb]{0,0,1}{#1}} $$

以下就是渲染方程，这个方程是基于辐射度量学的量化地描述渲染过程的公式。

$$ L_o(p, \omega) = L_e(p, \omega) + \int_{\Omega}f_r(\omega_o, \omega_i, p)L_i(p, \omega_o)cos {\theta}d\omega \\\\= L_e(p, \omega) + \int_{\Omega}f_r(\omega_o, \omega_i, p)L_i(p, \omega_o) (\omega_i \cdot n) d\omega $$

Cook-Torrance BRDF基于Microfacet模型得到了由Fresnel Formula, Normal Distribution Function 和 Geometry Term 组成的 PBR 的 BRDF.

它大概长这样

$$ f_r(w_o, w_i, p) = \frac{\red{F(\omega_o, h)}\blue{G(\omega_i, \omega_o)}\green{D(h)}}{4(n\cdot \omega_i)(n \cdot \omega_o)} $$

其中$h = normalize(\omega_i, \omega_o)$

Fresnel Formula

考虑到Fresnel Formula的复杂性，在工业上一般使用Fresnel-Schlick近似来模拟真正的Fresnel项。

$$ F_{Schlick}(\omega_o, h) = F_0 + (1 - F_0)(1 - (h \cdot \omega_o))^5 $$

其中$F_0$用来表示材质的基础反射率，也就是你沿法线方向观察物体，它的反射比重是怎样的。对于水面或塑料而言，你朝球面掠角看过去的话，Fresnel现象越明显，反光就越强。对于导体而言，几乎所有的颜色都由反射贡献，所以它们的$F_0$项差距很大。

那么我们用GLSL代码写出来是这样的

vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1. - F0) * pow(1. - cosTheta, 5.);
}

Normal Distribution Function

简称NDF，NDF做的就是类似之前我们在Blinn-Phong的Specular项里做的。但是这里引入了粗糙度这一物理量roughness.

之前工业界采用Beckmann的NDF，但是随着GGX的提出，长拖尾的NDF渲染起来更加自然，所以现代人们普遍使用GGX。

比较新的理论是有一种叫GTR（Generalized TR）的模型，它在GGX的基础上加入了新的参数$\gamma$ ，并且随着$\gamma$的变化可以和GGX与Beckmann等价。

下面是Trowbridge-Reitz GGX:

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

当粗糙度很低（也就是说表面很光滑）的时候，与中间向量取向一致的微平面会高度集中在一个很小的半径范围内。由于这种集中性，NDF最终会生成一个非常明亮的斑点。但是当表面比较粗糙的时候，微平面的取向方向会更加的随机。你将会发现与h向量取向一致的微平面分布在一个大得多的半径范围内，但是同时较低的集中性也会让我们的最终效果显得更加灰暗。

我们用GLSL去实现它

float GGX(vec3 N, vec3 H, float roughness)
{
    float seq_roughness = roughness * roughness;
    float NdotH = dot(H, N);
    float seq_NdotH = NdotH * NdotH;
    float nom = seq_roughness;
    float denom = (seq_NdotH * (seq_roughness - 1.) + 1.);
    denom = PI * denom * denom;
    return nom / denom;
}

Geometry Term

我们看下面两种情况，当一个平面凹凸不平（粗糙）或入射、反射光在掠角时，光有很明显的自遮挡的情况。这种情况我们需要一个项来修正。

另外一种理解方式是，当入射、反射光在掠角时，Cook-Torrance的分母将会无限接近于0，此时，我们的分子同样需要一个很小的项来中和它，来避免不自然的高光。

入射光被遮挡叫Shadowing，出射光被遮挡叫Masking.

另外，不同的法线分布函数对应的几何函数也互不相同，二者具有很强的依赖关系，需要搭配着一起使用。

与NDF类似，几何函数采用一个材料的粗糙度参数作为输入参数，粗糙度较高的表面其微平面间相互遮蔽的概率就越高。我们将要使用的几何函数是GGX与Schlick-Beckmann近似的结合体，因此又称为Schlick-GGX：

$$ G_{SchlickGGX}(\omega, n, k(\alpha)) = \frac{n \cdot \omega}{(n \cdot \omega)(1 - k(\alpha)) + k(\alpha)} $$

在这里$k(\alpha)$可以根据直接光照或IBL光照进行调整。

为了有效的估算Geometry Term，需要将观察方向（几何遮蔽(Geometry Obstruction)）和光线方向向量（几何阴影(Geometry Shadowing)）都考虑进去。我们可以使用史密斯法(Smith’s method)来把两者都纳入其中

$$ G(\omega_i, \omega_o, n, k(\alpha)) = G_{SchlickGGX}(\omega_i, n, k(\alpha)) \cdot G_{SchlickGGX}(\omega_o, n, k(\alpha)) $$

G项的值域是[0, 1]，这可以简单推出来（主要是1）

$$ G < \max_{u = \{i, o\} } \left[ \frac{n \cdot \omega_u}{(n \cdot \omega_u)(1 - k(\alpha)) + k(\alpha)} \right] ^2 \leq \frac{n \cdot \omega}{(n \cdot \omega)(1 - k(\alpha)) + (n \cdot \omega)k(\alpha)} = 1 $$

用GLSL实现一下，这里$G_{SchlickGGX}$可以视为只有两个参数，即

$$ G_{SchlickGGX}(\omega \cdot n, k(\alpha)) = G_{SchlickGGX}(\omega, n, k(\alpha)) $$

float SchlickGGX(float NdotV, float k)
{
    float nom = NdotV;
    float denom = NdotV * (1. - k) + k;
    return nom / denom;
}
float Smith(vec3 N, vec3 V, vec3 L, float k)
{
    float NdotV = max(0., dot(N, V));
    float NdotL = max(0., dot(N, L));
    return SchlickGGX(NdotV, k) * SchlickGGX(NdotL, k);
}

BRDF

考虑渲染方程，这里我们同样加入diffuse项

$$ L_o(p, \omega) = \int_{\Omega}\left( K_d\frac{c}{\pi} + K_s\frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}\right)L_i(p, \omega_o)cos {\theta} d\omega \\\\ = \int_{\Omega}\left( (1 - F)\frac{c}{\pi} + \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}\right)L_i(p, \omega_o)cos {\theta} d\omega $$

下面使用Cook-Torrance BRDF来渲染的样例。

out vec4 color;
in vec2 texCoord;
in vec3 normal;
in vec3 worldPos;

// uniform terms...
// function declaration

#define fr(name, begin, end) \
for(int name = begin; name < end; ++name)

#define max0f(x) max(0., x)

void main()
{
    vec3 N = Normal;
    vec3 V = normalize(camPos - worldPos);
    
    vec3 Lo = vec3(0.);
    fr(i, 0, LIGHTING_CNT)
    {
        vec3 L = normalize(lightPos[i] - worldPos);
        vec3 H = normalize(L + V);
        // attenuation
        float distance = length(lightPos[i] - worldPos);
        float attenuation = 1. / distance * distance;
        vec3 radiance = attenuation * lightColor[i];
        // brdf
        float brdf_n = GGX(N, H, roughness);
        float brdf_g = Smith(N, V, L, roughness);
        vec3 brdf_f = fresnelSchlick(max(dot(H, V)), F0);
        vec3 nom = brdf_n * brdf_f * brdf_g;
        float denom = 4. * max0f(dot(N, L)) * max0f(dot(N, V));
        vec3 specular = nom / denom;
        vec3 ks = F;
        vec3 kd = vec3(1.) - ks;
        kd *= 1.0 - metallic;
        float NdotL = max0f(dot(N, L));
        Lo += (kd * albedo / PI + specular) * radiance * NdotL;
    }
    // ...
}

这里我需要补充一下Attenuation和能量守恒

这个项其实是光线衰减项，想象一下烟花爆炸，理想情况下烟花是球形爆炸形状，也就是一个球壳上分布着相同的烟花粒子。因为亮度和烟花粒子呈线性关系，而球壳面积和半径的关系是平方关系，那么我们可以得到

$$ L \propto R^2 $$

相似的，我们可以得到Randiance衰减系数。

能量守恒即Ks和Kd项相加和必须为1. 因为导体大部分是反射项，我们让Ks项乘1减金属度。

IBL Environment Lighting

Diffuse Term

theory

我们对于环境光照的BRDF可以拆分成diffuse的和specular的。下面是diffuse的

$$ L_o(p, \omega) = \int_{\Omega}\left( K_d\frac{c}{\pi} \right)L_i(p, \omega_o)cos {\theta} d\omega \\\\ = K_d\frac{c}{\pi} \int_{\Omega}L_i(p, \omega_o)cos {\theta} d\omega \\\\ = K_d\frac{c}{\pi} \int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}} L_i(p, \phi, \theta) cos\theta sin\theta d\phi d\theta \\\\ \approx K_d\frac{c\pi}{n_1n_2} \sum_{\phi = 0}^{n_1}\sum_{\theta = 0}^{n_2} L_i(p, \phi, \theta) cos\theta sin\theta $$

在这里，我们就可以得到某个点的半球积分了。

pre-compute

我们如果对半球积分采样，就要遍历$\phi, \theta$

我们可以画一个cube，然后把上述积分转换成代码就好。

vec3 irradiance = vec3(0.);
vec3 up = vec3(0., 1., 0.);
vec3 right = normalize(cross(up, normal));
up = normalize(cross(normal, right));
float sampleDelta = .025;
float nrSamples = .0;
for (float phi = .0; phi < 2. * PI; phi += sampleDelta)
{
    for(float theta = .0; theta < .5 * PI; theta += sampleDelta)
    {
         vec3 tangentSample = vec3(sin(theta) * cos(phi),  sin(theta) * sin(phi), cos(theta));
         vec3 sampleVec = tangentSample.x * right + tangentSample.y * up + tangentSample.z * N;
         irradiance += texture(environmentCubeMap, sampleVec.rgb) * cos(theta) * sin(theta);
         nrSamples++;
    }
}
irradiance = PI * (irradiance / nrSamples);

那么我们就可以计算ambient项了，继续写之前的pbr-rendering.

uniform samplerCube irradianceMap;
// ...
void main()
{
    // pbr rendering.
    vec3 kS = fresnelSchlick(max(dot(N, V), 0.0), F0);
    vec3 kD = vec3(1.0) - kS;
    vec3 irradiance = texture(irradianceMap, N).rgb;
    vec3 diffuse    = irradiance * albedo;
    vec3 ambient    = (kD * diffuse) * ao; 
}

这里有些事很可疑，fresnel项需要一个中程向量H来计算相应的系数，但是考虑到环境光来自四面八方，没用一个单独的中程向量，但是为了模拟fresnel，我们用dot(N, V)模拟这一系数。

Fresnel Attenuation

thanks to Sébastien Lagarde.

这里基于一个理论，当斜视时，经过过滤和未经过滤的立方体贴图造成的光照差距是明显的，后者对多个方向的入射光进行平均并且用单一的Fresnel项计算反射权重，因为未考虑粗糙度，表面反射率总是相对较高。

所以对比较粗糙的表面进行矫正：

vec3 fresnelSchlickRoughness(float cosTheta, vec3 F0, float roughness)
{
    return F0 + (max(vec3(1.0 - roughness), F0) - F0) * pow(clamp(1.0 - cosTheta, 0.0, 1.0), 5.0);
}

然后之前的代码改为：

vec3 kS = fresnelSchlickRoughness(max(dot(N, V), 0.0), F0, roughness); 
vec3 kD = vec3(1.0) - kS;
vec3 irradiance = texture(irradianceMap, N).rgb;
vec3 diffuse    = irradiance * albedo;
vec3 ambient    = (kD * diffuse) * ao; 

然后接续之前的shader, 当前点的颜色就是

color = Lo + ambient;

Specular Term

theory

然后再看 Specular Term

$$ L_o(p, \omega) = \int_{\Omega} \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}L_i(p, \omega_o)cos {\theta} d\omega $$

pre compute