内容参考闫令琪课程《games202-高质量实时渲染及作业4》、“202作业4代码模板”、花桑博客
实现效果
原理说明
PBR材质分两大流派:
- 微表面模型:假设物体表面在微观层面是凹凸不平的
- 迪士尼材质:将物理属性转化成对艺术设计人员友好的参数-提高生产率,向外暴露了10个参数
微表面模型的BRDF有能量损失,并不是真正的能量守恒,粗糙度越大,几何项G值越小导致物体越暗,原因是G项只考虑光线一次弹射,忽略了光线的多次弹射
f ( ω i , ω o ) = F ( ω i , h ) G ( ω i , ω o , h ) D ( h ) 4 ( n ⋅ ω i ) ( n ⋅ ω o ) f(\omega_i, \omega_o)=\frac{\mathbf{F}(\omega_i,h)\mathbf{G}(\omega_i,\omega_o,h)\mathbf{D}(h)}{4(n\cdot \omega_i)(n\cdot \omega_o)} f(ωi,ωo)=4(n⋅ωi)(n⋅ωo)F(ωi,h)G(ωi,ωo,h)D(h)
Kulla-Conty近似的原理是基于原有的BRDF,再加上一个多次弹射的能量损失项,看下图,上面一张中粗糙度越大能量损失的越多模型整体偏黑,下图中补充了Kulla-Conty项,粗糙度变大镜面反射慢慢变成漫反射,颜色趋于灰白,但是整体亮度没有变暗
kulla-conty公式
kulla-conty近似推导的整体流程
- 不考虑光线吸收(只有镜面反射项)积分计算出BRDF基础项能量 E ( μ o ) E(\mu_o) E(μo)
- 设计一个kulla-conty补充项 f m s f_{ms} fms,使得其积分等于 1 − E ( μ o ) 1 - E(\mu_o) 1−E(μo)
- 考虑光线吸收项KaTeX parse error: \tag works only in display equations
白色部分波长的光线被物体吸收后,剩余光线再反射出来就得到物体的颜色
1)对BRDF进行积分,得到镜面反射项的总能量,便于处理设光照L项为1
E ( μ o ) = ∫ 0 2 π ∫ 0 1 f ( μ o , μ i , ϕ ) μ i d μ i d ϕ E(\mu_o)=\int_{0}^{2\pi}\int_{0}^1f(\mu_o,\mu_i,\phi)\mu_i d\mu_i d\phi E(μo)=∫02π∫01f(μo,μi,ϕ)μidμidϕ
其中 μ i = s i n θ i \mu_i=sin\theta_i μi=sinθi,
2)设计能量补充项 f m s f_{ms} fms,避免对每一种出射方向都做补充, E a v g E_{avg} Eavg是一种近似的考量,即不考虑方向了,每个方向都使用平均能量
f m s ( μ o , μ i ) = ( 1 − E ( μ o ) ) ( 1 − E ( μ i ) ) π ( 1 − E a v g ) f_{ms}(\mu_o,\mu_i)=\frac{(1-E(\mu_o))(1-E(\mu_i))}{\pi(1-E_{avg})} fms(μo,μi)=π(1−Eavg)(1−E(μo))(1−E(μi))
E a v g = 2 ∫ 0 1 E ( μ ) μ d μ E_{avg}=2\int_0^1 E(\mu)\mu d\mu Eavg=2∫01E(μ)μdμ
3)完整的kulla-conty项
f m s ( μ o , μ i ) = F a v g ⋅ E a v g 1 − F a v g ( 1 − E a v g ) ⋅ ( 1 − E ( μ o ) ) ( 1 − E ( μ i ) ) π ( 1 − E a v g ) f_{ms}(\mu_o,\mu_i)=\frac{F_{avg}\cdot E_{avg}}{1-F_{avg}(1-E_{avg})}\cdot\frac{(1-E(\mu_o))(1-E(\mu_i))}{\pi(1-E_{avg})} fms(μo,μi)=1−Favg(1−Eavg)Favg⋅Eavg⋅π(1−Eavg)(1−E(μo))(1−E(μi))
所以, f m s ( μ o , μ i ) + B R D F f_{ms}(\mu_o,\mu_i) + BRDF fms(μo,μi)+BRDF 可以得到新的BRDF项,不管粗糙度多大都能保持能量守恒,更符合真实的光照效果
重点 & 代码说明
预计算E(μ),c++工程中实现
将brdf项积分预计算好,生成图片,横轴对应粗糙度,纵轴对应cosΘ。这张图对应的一种粗糙度可变的材质
注意,实际代码的实现和课程里讲的细节有点不一样
课程里换元法,换掉了 c o s θ s i n θ cos\theta sin\theta cosθsinθ,但是代码里还是基于 f r ∗ c o s θ ∗ d w f_r * cos\theta * dw fr∗cosθ∗dw积分
作业中要求用两种方式实现brdf的预积分,一种普通的采样,一种是重要性采样(UE4里也是基于重要性采样)
1)普通的采样
代码中实现
//Emu_MC.cpp
Vec3f IntegrateBRDF(Vec3f V, float roughness, float NdotV) {
float A = 0.0;
float B = 0.0;
float C = 0.0;
const int sample_count = 1024;
Vec3f N = Vec3f(0.0, 0.0, 1.0);
samplePoints sampleList = squareToCosineHemisphere(sample_count);
for (int i = 0; i < sample_count; i++) {
// TODO: To calculate (fr * ni) / p_o here
// Edit Start
Vec3f L = normalize(sampleList.directions[i]);
float pdf = sampleList.PDFs[i];
Vec3f H = normalize(V + L);
float NdotL = std::max(dot(N, L), 0.0f);
float NDF = DistributionGGX(N, H, roughness);
float G = GeometrySmith(roughness, NdotV, NdotL);
float F = 1.0f;
float mu = NdotL;
float numerator = NDF * G * F;
float denominator = 4.0 * NdotV * NdotL;
A = B = C += numerator / denominator / pdf * mu;
// Edit End
}
return {A / sample_count, B / sample_count, C / sample_count};
}
预计算 E a v g E_{avg} Eavg
fms中E(μ)项已经预计算出来了,还缺 E a v g E_{avg} Eavg,也是预计算出来,不过 E a v g E_{avg} Eavg是一维的,放到图片里存储,只能让每一行的值相同
//Eavg_MC.cpp
Vec3f IntegrateEmu(Vec3f V, float roughness, float NdotV, Vec3f Ei) {
return Ei * NdotV * 2.0f;
}
注意前面E的图中,x轴是rough,y轴是 c o s θ cos\theta cosθ
2)重要性采样(重要性采样,可以参考前面的文章)
Hammersly函数-地差异序列,使样本更均匀
//Emu_IS.cpp
Vec2f Hammersley(uint32_t i, uint32_t N) { // 0-1
uint32_t bits = (i << 16u) | (i >> 16u);
bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
float rdi = float(bits) * 2.3283064365386963e-10;
return {float(i) / float(N), rdi};
}
重要的代码来了,重要性采样的GGX实现,将采样的(0,1)值变成3维的方向向量,这里课程上没仔细讲,learnOpengl中也是一笔带过,我查阅了不少资料才弄明白每一个细节
//Emu_IS.cpp
Vec3f ImportanceSampleGGX(Vec2f Xi, Vec3f N, float roughness) {
float a = roughness * roughness;
//TODO: in spherical space - Bonus 1
float theta = atan(a * sqrt(Xi.x) / sqrt(1.0f - Xi.x));
float phi = 2.0 * PI * Xi.y;
//TODO: from spherical space to cartesian space - Bonus 1
float sinTheta = sin(theta);
float consTheta = cos(theta);
Vec3f H = Vec3f(cos(phi) * sinTheta, sin(phi) * sinTheta, consTheta);
//TODO: tangent coordinates - Bonus 1
Vec3f up = abs(N.z) < 0.999 ? Vec3f(0.0, 0.0, 1.0) : Vec3f(1.0, 0.0, 0.0);
Vec3f tangent = normalize(cross(up, N));
Vec3f bitangent = cross(N, tangent);
//TODO: transform H to tangent space - Bonus 1
Vec3f sampleVec = tangent * H.x + bitangent * H.y + N * H.z;
return normalize(sampleVec);
}
最后同样也是积分,注意菲涅尔项为1,只考虑NDF和G项即可,不过没明白,积分公式的分母中4怎么不见了?后面有时间再好好研究下这个公式。
//Emu_IS.cpp
Vec3f IntegrateBRDF(Vec3f V, float roughness) {
const int sample_count = 1024;
Vec3f Emu(0.0f);
Vec3f N = Vec3f(0.0, 0.0, 1.0);
for (int i = 0; i < sample_count; i++) {
Vec2f Xi = Hammersley(i, sample_count);
Vec3f H = ImportanceSampleGGX(Xi, N, roughness);
Vec3f L = normalize(H * 2.0f * dot(V, H) - V);
float NoL = std::max(L.z, 0.0f);
float NoH = std::max(H.z, 0.0f);
float VoH = std::max(dot(V, H), 0.0f);
float NoV = std::max(dot(N, V), 0.0f);
// Edit Start
// TODO: To calculate (fr * ni) / p_o here - Bonus 1
float G = GeometrySmith(roughness, NoV, NoL);
float weight = VoH * G / (NoV * NoH);
Emu += Vec3f(1.0, 1.0, 1.0) * weight;
// Split Sum - Bonus 2
}
return Emu / sample_count;
}
同样也能生成一张图,这里就不贴了
回到前端js代码框架中,进行实时渲染
GGX_E_LUT.png和GGX_Eavg_LUT.png拷贝到实时渲染端下的assets/ball目录下,打开engine.js中的注释项
//engine.js
//..
loadGLTF(renderer, 'assets/ball/', 'ball', 'KullaContyMaterial', Sphere0Transform, metallic, 0.15);
//..
loadGLTF(renderer, 'assets/ball/', 'ball', 'PBRMaterial', Sphere5Transform, metallic, 0.15);
PBRFragment.glsl,KullaContyFragment.glsl两个Shader中补充微表面BRD的实现
//PBRFragment.glsl、KullaContyFragment.glsl
float DistributionGGX(vec3 N, vec3 H, float roughness)
{
// TODO: To calculate GGX NDF here
float a = roughness*roughness;
float a2 = a*a;
float NdotH = max(dot(N, H), 0.0);
float NdotH2 = NdotH*NdotH;
float nom = a2;
float denom = (NdotH2 * (a2 - 1.0) + 1.0);
denom = PI * denom * denom;
return nom / max(denom, 0.0001);
float GeometrySchlickGGX(float NdotV, float roughness)
{
// TODO: To calculate Smith G1 here
float a = roughness;
float k = (a * a) / 2.0;
float nom = NdotV;
float denom = NdotV * (1.0 - k) + k;
return nom / denom;
}
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness)
{
// TODO: To calculate Smith G here
float NdotV = max(dot(N, V), 0.0);
float NdotL = max(dot(N, L), 0.0);
float ggx2 = GeometrySchlickGGX(NdotV, roughness);
float ggx1 = GeometrySchlickGGX(NdotL, roughness);
return ggx1 * ggx2;
}
vec3 fresnelSchlick(vec3 F0, vec3 V, vec3 H)
{
// TODO: To calculate Schlick F here
return F0 + (1.0 - F0) * pow(clamp(1.0 - max(dot(H, V), 0.0), 0.0, 1.0), 5.0);
}
}
//KullaContyFragment.glsl
vec3 MultiScatterBRDF(float NdotL, float NdotV)
{
vec3 albedo = pow(texture2D(uAlbedoMap, vTextureCoord).rgb, vec3(2.2));
vec3 E_o = texture2D(uBRDFLut, vec2(NdotL, uRoughness)).xyz;
vec3 E_i = texture2D(uBRDFLut, vec2(NdotV, uRoughness)).xyz;
vec3 E_avg = texture2D(uEavgLut, vec2(0, uRoughness)).xyz;
// copper
vec3 edgetint = vec3(0.827, 0.792, 0.678);
vec3 F_avg = AverageFresnel(albedo, edgetint);
// TODO: To calculate fms and missing energy here
// TODO: To calculate fms and missing energy here
vec3 F_ms = (1.0 - E_o) * (1.0 - E_i) / (PI * (1.0 - E_avg));
vec3 F_add = F_avg * E_avg / (1.0 - F_avg * (1.0 - E_avg));
return F_add * F_ms;
}
注意KullaContyMaterial.js 中代码有一处错误, 需要把uEavgFLut改成uEavgLut