微表面模型属于材质方面，因此需要修改材质类Material.hpp。
本片文章补充了https://blog.csdn.net/qq_36242312/article/details/116307626?spm=1001.2014.3001.5502博客中的一些细节。
已知：
求:
基于物理的渲染（PBR）
f
s
=
D
∗
G
∗
F
4
∗
(
n
⃗
⋅
l
⃗
)
∗
(
n
⃗
⋅
v
⃗
)
f_{s}=\frac{D * G * F}{4 *(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})}
fs=4∗(n
⋅l
)∗(n
⋅v
)D∗G∗F
微平面的法线分布函数D(m)描述了微观表面上的表面法线m的统计分布。
业界较为主流的法线分布函数是GGX（Trowbridge-Reitz），因为具有更好的高光长尾。
D
(
n
,
h
,
α
)
=
α
2
π
(
(
n
⋅
h
)
2
(
α
2
−
1
)
+
1
)
2
D(n, h, \alpha) = \frac{\alpha^{2}}{\pi \left((n \cdot h)^{2}\left(\alpha^{2}-1\right)+1\right)^{2}}
D(n,h,α)=π((n⋅h)2(α2−1)+1)2α2
当
α
=
1
\alpha=1
α=1时，
D
=
1
π
D = \frac 1 \pi
D=π1,光线向半球面均匀发散;
当
α
=
0
\alpha=0
α=0时，只有当
n
=
h
时
，
D
≠
0
n=h时，D \ne0
n=h时，D=0 。
如果想要理解函数是怎么来的，建议查阅相关论文，作为小白的我，已放弃挣扎~~
float D = DistributionGGX(N, H, roughness);
根据公式计算D值：
a为
α
\alpha
α
NdotH为
N
和
H
N和H
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
}
目前较为常用的是其中最为简单的形式，分离遮蔽阴影（Separable Masking and Shadowing Function）。
该形式将几何项G分为两个独立的部分：光线方向（light） 和 视线方向（view），并对两者用相同的分布函数来描述。
其中UE4的方案是Schlick-GGX，即基于Schlick近似，将k映射为 k = α 2 k = \frac {\alpha}{2} k=2α ,去匹配GGX Smith方程：
α
=
(
r
o
u
g
h
t
n
e
s
s
+
1
2
)
2
−
−
−
−
−
−
−
−
−
−
−
−
k
=
α
2
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
G
d
i
r
e
c
t
i
o
n
(
n
,
v
,
k
)
=
n
⋅
v
(
n
⋅
v
)
(
1
−
k
)
+
k
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
G
(
n
,
v
,
l
,
k
)
=
G
v
i
e
w
(
n
,
v
,
k
)
G
l
i
g
h
t
(
n
,
l
,
k
)
\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)
α=(2roughtness+1)2−−−−−−−−−−−−k=2α−−−−−−−−−−−−−−−−Gdirection(n,v,k)=(n⋅v)(1−k)+kn⋅v−−−−−−−−−−−−−−−−−−−−−G(n,v,l,k)=Gview(n,v,k)Glight(n,l,k)
G
(
n
,
v
,
l
,
k
)
G(n, v, l, k)
G(n,v,l,k)返回遮蔽阴影函数的计算结果。
因此 G ( n , v , l , k ) ∈ [ ( n ⋅ v ) ( n ⋅ l ) , 1 ] G(n, v, l, k) \in [(n \cdot v)(n \cdot l),1] G(n,v,l,k)∈[(n⋅v)(n⋅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）项，业界方案一般都采用Schlick的Fresnel近似，因为计算成本低廉，而且精度足够：
F
s
c
h
l
i
c
k
(
v
,
h
)
=
F
0
+
(
1
−
F
0
)
(
1
−
(
v
⋅
h
)
)
5
F_{schlick}(v,h) = F_0 + (1-F_0)(1-(v \cdot h))^5
Fschlick(v,h)=F0+(1−F0)(1−(v⋅h))5
但是在项目中没有使用近似计算，而是使用了Fresnel方程。
分别对应的菲涅尔方程为：
详细内容看菲涅尔方程（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;
}
Cook-Torrance BRDF兼有漫反射和镜面反射两个部分：
f r = k d ∗ f d + k s ∗ f s f_r=k_d*f_{d}+k_s*f_{s} fr=kd∗fd+ks∗fs
f s = D ∗ G ∗ F 4 ∗ ( n ⃗ ⋅ l ⃗ ) ∗ ( n ⃗ ⋅ v ⃗ ) f_{s}=\frac{D * G * F}{4 *(\vec{n} \cdot \vec{l}) *(\vec{n} \cdot \vec{v})} fs=4∗(n ⋅l )∗(n ⋅v )D∗G∗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);
specular 即为
f
s
f_s
fs 。
f d = c π f_d = \frac c \pi fd=πc
此处我也没搞懂，以后再填坑 。
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;
}
}
}
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;
}
}
}
在Intersection getIntersection(Ray ray)函数中修改圆的相交判定，因为圆的相交判定精度不够。
若不修改，程序在运行中会报错，因为当
t
0
=
0
t0=0
t0=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) ,则代码不会在运行中报错，但是会生成有黑点的图像。
若设置精度为t>0.1,生成的图像为：
可以看出t>0.1和t>0.5的差别并不大。那为什么（t!=0或t>0)这样的条件会导致图像黑点呢？？？
下面我们来分析原因：
光线与球体求交公式使用了一般情况下的求交公式（在该项目中可以优化求交公式。因为ray.origin即光源位置；ray.direction为单位向量）。
而求交公式获得的交点有偏差，如下图：
实际光线与球面相交返回的交点在球面附近。
根据程序流程，找到某个物体的交点后，要将该交点和某个光源连线，判断该点光源是否可见。
因此会出现上图中的有黑点的球体。
在交点和光源的连线和球面相交时，若球面距离弹射点过近，则认为球面和光线不相交。这样便可避免上方的错误。
这便是在Intersection getIntersection(Ray ray)函数中要限定t>某个值的意义！！！
在下面函数中添加返回值，不然会报错。
Vector3f evalDiffuseColor(const Vector2f &st)const {
//return m->getColor();
return {};
}
最后祝大家新春快乐
GAMES101-现代计算机图形学学习笔记（作业07）
从零开始学图形学：写一个光线追踪渲染器（二）——微表面模型与代码实现
LearnOpenGL_CN 理论
【基于物理的渲染（PBR）白皮书】（一） 开篇：PBR核心知识体系总结与概览
【基于物理的渲染（PBR）白皮书】（四）法线分布函数相关总结
【基于物理的渲染（PBR）白皮书】（五）几何函数相关总结
计算机图形学——学习记录三（光线相交）
菲涅尔方程（Fresnel Equation）