Games101,作业7（微表面模型）

Elsa的迷弟

已于 2022-03-30 12:56:40 修改

阅读量2.3k

点赞数 9

分类专栏： Games101（图形学基础）文章标签： games101 c++ 微表面模型

于 2022-01-27 15:11:25 首次发布

本文链接：https://blog.csdn.net/weixin_44518102/article/details/122698851

版权

Games101（图形学基础）专栏收录该内容

21 篇文章 10 订阅

订阅专栏

微表面模型

微表面模型属于材质方面，因此需要修改材质类Material.hpp。

本片文章补充了https://blog.csdn.net/qq_36242312/article/details/116307626?spm=1001.2014.3001.5502博客中的一些细节。

BRFD

已知：

入射方向 $w i$
出射方向 $w o$
法线 $N$

求:

光线反射的（RGB）反射率 $f_r$ 。

微表面模型的BRFD

基于物理的渲染（PBR）
$f_{s}=\frac{D * G * F}{4 *(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})}$

D 法线分布函数：某个点上，法线指向指定方向上的概率。
G 几何函数：光线由于遮蔽、阴影，所以部分光线被阻挡；G即为光线没有被阻挡的概率。
F Fresnel方程：表示一束光线打在某个方向的平面上的反射的能量比率。
$*(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})}$ 为矫正量

微表面模型需要用到的参数:

$V$ ：光线的反射方向（单位向量）
$I$ ：光线的入射方向的反方向（单位向量）
$m = h$ ： $V$ 和 $I$ 的中间向量（单位向量）
$N$ ：为宏观平面法线（单位向量）
$\alpha$ ：为粗糙度 $\in[0,1]$ 。

法线分布函数 D

微平面的法线分布函数D(m)描述了微观表面上的表面法线m的统计分布。

业界较为主流的法线分布函数是GGX（Trowbridge-Reitz），因为具有更好的高光长尾。
$\alpha) = \frac{\alpha^{2}}{\pi \left((n \cdot h)^{2}\left(\alpha^{2}-1\right)+1\right)^{2}}$

$\alpha$ 为粗糙度 $\in[0,1]$ 。
$n$ 为宏观平面法线
$h$ 为微观平面法线，即 $V$ 和 $I$ 的中间向量。

当 $\alpha=1$ 时， $\frac 1 \pi$ ,光线向半球面均匀发散;
当 $\alpha=0$ 时，只有当 $\ne0$ 。

如果想要理解函数是怎么来的，建议查阅相关论文，作为小白的我，已放弃挣扎~~

代码：

float D = DistributionGGX(N, H, roughness);

根据公式计算D值：

a为 $\alpha$
NdotH为 $N 和 H$ 的点乘，必须大于0。
nom为分子
denom为分母（分母不能为0）。

 float DistributionGGX(Vector3f N, Vector3f H, float roughness)
    {
        float a = roughness * roughness;
        float a2 = a * a;
        float NdotH = std::max(dotProduct(N, H), 0.0f);
        float NdotH2 = NdotH * NdotH;

        float nom = a2;
        float denom = (NdotH2 * (a2 - 1.0) + 1.0);
        denom = M_PI * denom * denom;

        return nom / std::max(denom, 0.0000001f); // prevent divide by zero for roughness=0.0 and NdotH=1.0
    }

几何函数 G

目前较为常用的是其中最为简单的形式，分离遮蔽阴影（Separable Masking and Shadowing Function）。

该形式将几何项G分为两个独立的部分：光线方向（light） 和 视线方向（view），并对两者用相同的分布函数来描述。

其中UE4的方案是Schlick-GGX，即基于Schlick近似，将k映射为 $\frac {\alpha}{2}$ ,去匹配GGX Smith方程：

$\alpha = (\frac {roughtness+1}{2})^2 \\------------\\ k = \frac {\alpha}{2} \\----------------\\ G_{direction}(n, v, k)=\frac{n \cdot v}{(n \cdot v)(1-k)+k} \\---------------------\\ G(n, v, l, k)=G_{view}(n, v, k) G_{light}(n, l, k)$
$G (n, v, l, k)$ 返回遮蔽阴影函数的计算结果。

粗糙度 $\alpha = 0$ 时,完全光滑，光线通过率为1.
粗糙度 $\alpha = 1$ 时，最粗糙，光线通过率为 $\cdot v$

因此 $\in [(n \cdot v)(n \cdot l),1]$

代码：

float G = GeometrySmith(N, V, L, roughness);

N为宏观表面法线
V为光线反射方向（可理解为视口的方向）
L为光线进入入的反方向（可理解为光源的方向）
roughness为粗糙度。

    float GeometrySchlickGGX(float NdotV, float k)
    {
        float nom = NdotV;
        float denom = NdotV * (1.0 - k) + k;

        return nom / denom;
    }

    float GeometrySmith(Vector3f N, Vector3f V, Vector3f L, float roughness)
    {
        float r = (roughness + 1.0);
        float k = (r * r) / 8.0;
        float NdotV = std::max(dotProduct(N, V), 0.0f);
        float NdotL = std::max(dotProduct(N, L), 0.0f);
        float ggx2 = GeometrySchlickGGX(NdotV, k);
        float ggx1 = GeometrySchlickGGX(NdotL, k);

        return ggx1 * ggx2;
    }

Fresnel方程 F

对于菲涅尔（Fresnel）项，业界方案一般都采用Schlick的Fresnel近似，因为计算成本低廉，而且精度足够：
$F_{schlick}(v,h) = F_0 + (1-F_0)(1-(v \cdot h))^5$

但是在项目中没有使用近似计算，而是使用了Fresnel方程。

$R_S$ 为S偏振光
$R_P$ 为P偏振光

分别对应的菲涅尔方程为：
在这里插入图片描述

详细内容看菲涅尔方程（Fresnel Equation）

代码

i为入射角度
t为折射角度
ior为平面物体的折射率（法线反向空间中介质的折射率）
kr为需要计算的反射率（即反射的光线所占的比率）。

在这里插入图片描述

// 计算 fresnel 系数: F
float F;
float etat = 1.85;
fresnel(wi, N, etat, F);

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

Cook-Torrance BRDF兼有漫反射和镜面反射两个部分：

$f_r=k_d*f_{d}+k_s*f_{s}$

微表面镜面反射项

$f_{s}=\frac{D * G * F}{4 *(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})}$

Vector3f nominator = D * G * F;
float denominator = 4 * std::max(dotProduct(N, V), 0.0f) \
 					  * std::max(dotProduct(N, L), 0.0f);
Vector3f specular = nominator / std::max(denominator, 0.001f);

specular 即为 $f_s$ 。

微表面漫反射项

$f_d = \frac c \pi$

c表示表面颜色
除以π是为了对漫反射光进行标准化，因为前面含有BRDF的积分方程是受π影响的。

此处我也没搞懂，~~以后再填坑~~ 。

其他项

Kd：为折射入物体内的光线从内部弹射出来后，没有被损失的能量以漫反射的方式分散到各个方向。
Ks：为反射的颜色分量，与物体表面颜色有关。

代码

// 能量守恒
float ks_ = F;//反射比率
float kd_ = 1.0f - ks_;//折射比率

Vector3f diffuse = 1.0f / M_PI;//漫反射系数

// 因为在 specular 项里已经考虑了反射部分的比例：F。所以反射部分不需要再乘以 ks_ 
return Ks * specular + kd_ * Kd * diffuse;

完整代码

private:
    float DistributionGGX(Vector3f N, Vector3f H, float roughness)
    {
        float a = roughness * roughness;
        float a2 = a * a;
        float NdotH = std::max(dotProduct(N, H), 0.0f);
        float NdotH2 = NdotH * NdotH;

        float nom = a2;
        float denom = (NdotH2 * (a2 - 1.0) + 1.0);
        denom = M_PI * denom * denom;

        return nom / std::max(denom, 0.0000001f); // prevent divide by zero for roughness=0.0 and NdotH=1.0
    }

    float GeometrySchlickGGX(float NdotV, float k)
    {
        float nom = NdotV;
        float denom = NdotV * (1.0 - k) + k;
        return nom / denom;
    }

    float GeometrySmith(Vector3f N, Vector3f V, Vector3f L, float roughness)
    {
        float r = (roughness + 1.0);
        float k = (r * r) / 8.0;
        float NdotV = std::max(dotProduct(N, V), 0.0f);
        float NdotL = std::max(dotProduct(N, L), 0.0f);
        float ggx2 = GeometrySchlickGGX(NdotV, k);
        float ggx1 = GeometrySchlickGGX(NdotL, k);

        return ggx1 * ggx2;
    }

Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        {
            // calculate the contribution of diffuse   model
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) {
                Vector3f diffuse = Kd / M_PI;
                return diffuse;
            }
            else
                return Vector3f(0.0f);
            break;
        }
        case Microfacet://微表面材质的BRDF
        {
            // Disney PBR 方案
            float cosalpha = dotProduct(N, wo);
            if (cosalpha > 0.0f) {
                float roughness = 0.40;

                Vector3f V = -wi;
                Vector3f L = wo;
                Vector3f H = normalize(V + L);

                // 计算 distribution of normals: D
                float D = DistributionGGX(N, H, roughness);

                // 计算 shadowing masking term: G
                float G = GeometrySmith(N, V, L, roughness);

                // 计算 fresnel 系数: F
                float F;
                float etat = 1.85;
                fresnel(wi, N, etat, F);

                Vector3f nominator = D * G * F;
                float denominator = 4 * std::max(dotProduct(N, V), 0.0f) * std::max(dotProduct(N, L), 0.0f);
                Vector3f specular = nominator / std::max(denominator, 0.001f);

                // 能量守恒
                float ks_ = F;//反射比率
                float kd_ = 1.0f - ks_;//折射比率

                Vector3f diffuse = 1.0f / M_PI;

                // 因为在 specular 项里已经考虑了反射部分的比例：F。所以反射部分不需要再乘以 ks_ 
                //Ks为镜面反射项，Kd为漫反射项。
                return Ks * specular + kd_ * Kd * diffuse;
            }
            else
                return Vector3f(0.0f);
            break;
        }
    }
}

其他需要修改的函数

in Material.hpp

Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
    switch(m_type){
    case DIFFUSE:
    case Microfacet:
        {
            // uniform sample on the hemisphere在半球上均匀采样
            float x_1 = get_random_float(), x_2 = get_random_float();
            //z∈[0,1]，是随机半球方向的z轴向量
            float z = std::fabs(1.0f - 2.0f * x_1);
            //r是半球半径随机向量以法线为旋转轴的半径
            //phi是r沿法线旋转轴的旋转角度
            float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;//phi∈[0,2*pi]
            Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);//半球面上随机的光线的弹射方向
            return toWorld(localRay, N);//转换到世界坐标
            
            break;
        }
    }
}

float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
    switch(m_type){
        case DIFFUSE:
        case Microfacet:
        {
            // uniform sample probability 1 / (2 * PI)
            if (dotProduct(wo, N) > 0.0f)
                return 0.5f / M_PI;
            else
                return 0.0f;
            break;
        }
    }
}

in sphere.hpp

1.计算误差导致的交点偏离问题

在Intersection getIntersection(Ray ray)函数中修改圆的相交判定，因为圆的相交判定精度不够。
若不修改，程序在运行中会报错，因为当 $t 0 = 0$ 时，有可能会造成除0错误。

	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;
    }

生成的图片为：
在这里插入图片描述

若设置为 if (t0 != 0.0) ,则代码不会在运行中报错，但是会生成有黑点的图像。
blackSphere
若设置精度为t>0.1,生成的图像为：
在这里插入图片描述

可以看出t>0.1和t>0.5的差别并不大。那为什么（t!=0或t>0)这样的条件会导致图像黑点呢？？？

下面我们来分析原因：

光线与球体求交

光线与球体求交公式使用了一般情况下的求交公式（在该项目中可以优化求交公式。因为ray.origin即光源位置；ray.direction为单位向量）。

而求交公式获得的交点有偏差，如下图：请添加图片描述
实际光线与球面相交返回的交点在球面附近。
根据程序流程，找到某个物体的交点后，要将该交点和某个光源连线，判断该点光源是否可见。

若交点在球内部，则若按照一般情况判断，交点和光源的连线必然会被该球体遮蔽，光路可能在球体内部弹射（这不符合实际）。
若交点在球外部，则该点可以被光源或其他弹射光线照亮。

因此会出现上图中的有黑点的球体。

算法纠正：

在交点和光源的连线和球面相交时，若球面距离弹射点过近，则认为球面和光线不相交。这样便可避免上方的错误。

这便是在Intersection getIntersection(Ray ray)函数中要限定t>某个值的意义！！！

2.无返回值

在下面函数中添加返回值，不然会报错。

Vector3f evalDiffuseColor(const Vector2f &st)const {
        //return m->getColor();
        return {};
}

最终渲染效果

请添加图片描述

最后祝大家新春快乐
请添加图片描述

参考资料：

GAMES101-现代计算机图形学学习笔记（作业07）
从零开始学图形学：写一个光线追踪渲染器（二）——微表面模型与代码实现
 LearnOpenGL_CN 理论
 【基于物理的渲染（PBR）白皮书】（一）开篇：PBR核心知识体系总结与概览
 【基于物理的渲染（PBR）白皮书】（四）法线分布函数相关总结
 【基于物理的渲染（PBR）白皮书】（五）几何函数相关总结
 计算机图形学——学习记录三（光线相交）
菲涅尔方程（Fresnel Equation）

Elsa的迷弟

关注

9
点赞
踩
13

收藏

觉得还不错? 一键收藏
打赏
7
评论
Games101,作业7（微表面模型）

微表面模型微表面模型属于材质方面，因此需要修改材质类Material.hpp。本片文章补充了https://blog.csdn.net/qq_36242312/article/details/116307626?spm=1001.2014.3001.5502博客中的一些细节。BRFD已知：入射方向wiwiwi出射方向wowowo法线NNN求:光线反射的（RGB）反射率frf_rfr。微表面模型的BRFD基于物理的渲染（PBR）fs=D∗G∗F4∗(n⃗⋅l⃗)∗(n⃗⋅v⃗
复制链接

扫一扫