微表面模型 | Portfolio & Blog

参考资料：

https://learnopengl.com/PBR/Theory

https://zhuanlan.zhihu.com/p/606074595

本文的内容包括微表面模型的介绍以及代码实现，篇幅较长，如果不想看原理可以直接移步到代码实现。

微表面模型的定义

微表面模型，也可以认为是微表面材质，被广泛应用到所有的基于物理渲染技术（Physically Based Rendering，PBR）。想要实现PBR就绕不开这个材质。顺带一提，实现PBR还需要满足以下三个条件：

基于微表面模型
能量守恒原则
使用基于物理的BRDF

需要注意的是，PBR仍然是对现实光照和材质的模拟，只是他用到了基于物理的理论，而Blinn-Phong模型不是基于物理的，因为他缺少了能量守恒的条件，并且其用于描述光线与物体表面的作用使用到的数学模型也是基于经验的。

应用了微表面材质的物体，在远距离观察下表面是平的，凑近了看则是有细节的，表现为带有颗粒状以及镜面反射。这是因为物体表面被分为了一个个微小的表面，每个微表面具有独立的法线，能够以不规律的方式反射光线。

如果微表面的法线较为集中，则表现为Glossy一样的材质，如果微表面的法线比较分散，则表现为Diffuse，漫反射材质

BRDF

BRDF，可以被认为是用来描述材质的物理模型，不同的材质具有不同的物理表达式。

本文实现部分使用到的BRDF是Cook-Torrance BRDF，它包括漫反射和镜面反射两部分，如式（1）所示：
$$
f_r = k_d f_{lambert}+k_s f_{cook-torrance} \tag{1}
$$
其中Kd是入射光能量的折射比例，Ks 是反射比例，Kd和Ks之间满足Kd = 1 - Ks的关系，而Ks其实就是 DFG 中的F，也就是菲涅尔项（在Cook-Torrance 高光反射的部分会涉及）

这是因为反射光和折射光之间是互斥的，反射的任何光能量都不再被材料本身吸收。也就是说光线在物体表面发生了反射带走了一部分能量，而剩下的进入表面的就是入射光的能量，两者加起来不会超过1，满足了能量守恒。

Lambertian 漫反射

BRDF中的漫反射部分使用的是Lambertian漫反射，如式（2）所示：
$$
f_{lambert} = \frac{c}{\pi} \tag{2}
$$
其中c是反照率（albedo）或表面颜色（漫反射表面纹理）。除以 π 的目的是为了将漫反射光归一化。

Cook-Torrance 高光反射

Cook-Torrance 高光反射的表达式（3）如下：
$$
f(i,o) = \frac{F(i,h)G(i,o,h)D(h)}{4(n,i)(n,o)} \tag{3}
$$
其中，F是菲涅尔项，G是阴影遮罩项，D是法线分布项，以下是他们具体描述：

D (normal Distribution function，NDF) - 法线分布函数：用来近似微平面对半程向量的对准程度，受表面粗糙度影响；这是近似微平面的主要函数。
F (Fresnel equation) - 菲涅尔方程：菲涅尔方程描述了不同表面角度处的表面反射比。
G (Geometry function) - 几何函数：描述微平面的自遮挡特性。当表面相对粗糙时，表面的微平面可以遮挡其他微平面，减少表面反射的光线。

下图描述了这个方程中每个物理量的含义，h是半程向量，它是视线方向和光线方向一半方向上的单位向量。而在使用光线追踪的渲染中，wi也可以表示为入射光，wo表示为出射光。

需要注意的是，所有的方向都是单位向量。

法线分布函数（NDF）

法线分布函数有许多种，本文使用的是Trowbridge-Reitz GGX的一种：
$$
NDF_{GGXTR}(n, h, \alpha) = \frac{\alpha^2}{\pi ((n \cdot h)^2 (\alpha^2 - 1) + 1)^2} \tag{4}
$$
其中h是半程向量，α是表面粗糙度。

当粗糙度较低，表面较光滑时，大量微平面都集中对齐到一个小半径的半程向量。由于这种高度集中，物体会显示出一个非常明亮的点，就像下图这个球形材质的外观展示出来的一样，下图展示的是粗糙度为0.1的结果。

在表面粗糙度较高，也就是表面比较粗糙的情况下，微表面的法线分布更加分散，结果也更像漫反射的结果，下图展示的是粗糙度为0.8的结果。

几何函数

与NDF类似，几何函数将材质的粗糙度参数作为输入，表面越粗糙，微表面细节相互遮挡的概率就越高。

本文使用的几何函数是GGX和Schlick-Beckmann逼近的组合，称为Schlick-GGX，如式（5）所示：
$$
G_{ShlickGGX}(n,v,k) = \frac{n \cdot v}{(n \cdot v)(1-k)+k} \tag{5}
$$

其中，k使用的是和粗糙度有关的基于直接光照的系数，还有一种系数用于IBL照明，有兴趣的可以自行了解。
$$
K_{direct} = \frac{(\alpha + 1)^2}{8} \ \tag{6}
$$

需要同时考虑视线方向（几何遮挡）和光照方向向量（几何阴影），可以使用史密斯方法同时考虑这两个因素。具体实现为使用Schlick-GGX （5）作为 G_sub 得到视线方向和光线方向的几何函数，将它相乘即可得到外面最后的表达式，如式（7）所示：
$$
G(n,v,l,k) = G_{sub}(n,v,k)G_{sub}(n,l,k) \tag{7}
$$

菲涅尔方程

菲涅尔方程描述了在我们观察表面时反射光与折射光的比例，这个比例随着我们观察表面的角度而变化。

由于菲涅尔方程比较复杂，实际应用都是使用近似方程菲涅尔-施里克（Fresnel-Schlick），如式（8）所示：
$$
F_{Schlick}(h,v,F_0) = F_0 + (1-F_0)(1-(h \cdot v))^5 \tag{8}
$$

其中，F0 表示表面的基础反射率，我们使用反射率IOR来计算，可以根据你需要渲染的材质寻找相应的IOR。

C++代码实现

Lambertian 漫反射

由于Kd = 1 - Ks的关系，而Ks其实就是 DFG 中的F，因此代码如下：

1	Vector3f diffuse = (Vector3f(1.0f) - F) * color / M_PI;

Cook-Torrance 高光反射

NDF

根据公式敲就可以了，需要注意的是，n · h需要在0和他们之间取最大值，因为小于0是没有意义的，此外，需要对分母进行限制，防止过大的值。

代码如下：

float DistributionGGX(Vector3f N, Vector3f H, float Roughness)
{
    float alpha2 = Roughness * Roughness;
    float NdotH = dotProduct(N, H) > 0.0f ? dotProduct(N, H) : 0.0f;
    float NdotH2 = NdotH * NdotH;

    float nom = alpha2;
    float denom = (NdotH2 * (alpha2 - 1.0) + 1.0);
    denom = M_PI * denom * denom;
    if (denom < 0.001)
        return 1.0f;
    else
        return nom / denom;
}

几何函数

由于使用的是直接光照，k也要使用直接光照的系数。然后根据公式敲即可，代码如下：

float GeometrySmith(Vector3f N, Vector3f I, Vector3f L, float Roughness)
{
	float k = (Roughness + 1) * (Roughness + 1) / 8;

	float NdotV = dotProduct(N, I) > 0.0f ? dotProduct(N, I) : 0.0f;
	float NdotL = dotProduct(N, L) > 0.0f ? dotProduct(N, L) : 0.0f;

	float ggx1 = NdotV / (NdotV * (1.0 - k) + k);
	float ggx2 = NdotL / (NdotL * (1.0 - k) + k);

	return ggx1 * ggx2;
}

菲涅尔项

这里使用了作业框架中的实现的代码，作业框架的实现是直接使用准确的方程，如下图所示：

具体代码如下：

void fresnel(const Vector3f &I, const Vector3f &N, const float &ior, float &kr) const
{
    float cosi = clamp(-1, 1, dotProduct(I, N));
    float etai = 1, etat = ior;
    if (cosi > 0) {  std::swap(etai, etat); }
    // Compute sini using Snell's law
    float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));
    // Total internal reflection
    if (sint >= 1) {
        kr = 1;
    }
    else {
        float cost = sqrtf(std::max(0.f, 1 - sint * sint));
        cosi = fabsf(cosi);
        float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));
        float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));
        kr = (Rs * Rs + Rp * Rp) / 2;
    }
    // As a consequence of the conservation of energy, transmittance is given by:
    // kt = 1 - kr;
}

BRDF

最后将上面的代码整合在一起，得到了BRDF的代码，这里需要注意wi的方向，wi取正数表示的是入射光线的方向，而wi取负表示视线方向。

case MICROFACET:
{
    float cosalpha = dotProduct(N, wo);
    if (cosalpha > 0.0f) {

        Vector3f H = (-wi + wo).normalized();
        float Roughness = 0.25f;

        float D = DistributionGGX(N, H, Roughness);
        float G = GeometrySmith(N, -wi, wo, Roughness);
        float F;
        fresnel(wi, N, ior, F);

        Vector3f diffuse = (Vector3f(1.0f) - F) * Kd / M_PI;
        Vector3f specular;
        float divisor = 4 * dotProduct(wo, N) * dotProduct(-wi, N);
        //防止分母趋于0的情况导致值过大
        if (divisor < 0.001)
        {
            specular = Vector3f(1);
        }
        else
        {
            //specular = Ks * (D * F * G) / divisor;
            specular = (D * F * G) / divisor;
        }

        return diffuse + specular;
    }
    else
        return Vector3f(0.0f);
    break;
}

常见问题

如果渲染物体出现了黑色噪点过多的情况，有可能是判断光线求交中的精度问题导致没有得到正确的颜色值，可以修改对应的求交判断阈值，比如我这里赋予球体一个微表面材质，但是渲染结果中球体的正面出现了黑色的噪点，需要去Sphere.hpp中的判断相交的函数中修改判断求交的阈值精度。

Intersection getIntersection(Ray ray){
    Intersection result;
    result.happened = false;
    Vector3f L = ray.origin - center;
    float a = dotProduct(ray.direction, ray.direction);
    float b = 2 * dotProduct(ray.direction, L);
    float c = dotProduct(L, L) - radius2;
    float t0, t1;
    if (!solveQuadratic(a, b, c, t0, t1)) return result;
    if (t0 < 0) t0 = t1;
    if (t0 < 0) return result;
    //黑色噪点太多，修改这里的判断阈值
    if (t0 > 0.5)
    {
        result.happened = true;

        result.coords = Vector3f(ray.origin + ray.direction * t0);
        result.normal = normalize(Vector3f(result.coords - center));
        result.m = this->m;
        result.obj = this;
        result.distance = t0;
    }

    return result;

}

渲染结果

以下渲染结果采用的spp为1024，从上到下，从左到右的粗糙度依次为0.0、0.25、0.5以及0.75。

作业7-微表面材质