GAMES202实时渲染(4)-Kulla-Conty BRDF

内容参考闫令琪课程《games202-高质量实时渲染及作业4》、“202作业4代码模板”、花桑博客

实现效果

原理说明

PBR材质分两大流派：

• 微表面模型：假设物体表面在微观层面是凹凸不平的

• 迪士尼材质：将物理属性转化成对艺术设计人员友好的参数-提高生产率，向外暴露了10个参数

微表面模型的BRDF有能量损失，并不是真正的能量守恒，粗糙度越大，几何项G值越小导致物体越暗，原因是G项只考虑光线一次弹射，忽略了光线的多次弹射
$$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)}$$

Kulla-Conty近似的原理是基于原有的BRDF，再加上一个多次弹射的能量损失项，看下图，上面一张中粗糙度越大能量损失的越多模型整体偏黑，下图中补充了Kulla-Conty项，粗糙度变大镜面反射慢慢变成漫反射，颜色趋于灰白，但是整体亮度没有变暗

kulla-conty公式

kulla-conty近似推导的整体流程

• 不考虑光线吸收(只有镜面反射项)积分计算出BRDF基础项能量$E({\mu }_o)$
• 设计一个kulla-conty补充项$f_{ms}$，使得其积分等于$1-E({\mu }_o)$
• 考虑光线吸收项KaTeX parse error: \tag works only in display equations

白色部分波长的光线被物体吸收后，剩余光线再反射出来就得到物体的颜色

1）对BRDF进行积分，得到镜面反射项的总能量，便于处理设光照L项为1
$$E({\mu }_o)={\int }_0^{2\pi }{\int }_0^{1}f({\mu }_o,{\mu }_i,\varphi ){\mu }_{i}d{\mu }_{i}d\varphi$$
其中${\mu }_i=sin{\theta }_i$,

2）设计能量补充项$f_{ms}$，避免对每一种出射方向都做补充，$E_{avg}$是一种近似的考量，即不考虑方向了，每个方向都使用平均能量
$$f_{ms}({\mu }_o,{\mu }_i)=\frac{(1-E({\mu }_o))(1-E({\mu }_i))}{\pi (1-E_{avg})}$$

$$E_{avg}=2{\int }_0^{1}E(\mu )\mu d\mu$$

3）完整的kulla-conty项
$$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})}$$

所以，$f_{ms}({\mu }_o,{\mu }_i)+BRDF$ 可以得到新的BRDF项，不管粗糙度多大都能保持能量守恒，更符合真实的光照效果

重点 & 代码说明

预计算E(μ)，c++工程中实现

将brdf项积分预计算好，生成图片，横轴对应粗糙度，纵轴对应cosΘ。这张图对应的一种粗糙度可变的材质

注意，实际代码的实现和课程里讲的细节有点不一样

课程里换元法，换掉了$cos\theta sin\theta$，但是代码里还是基于$f_r\ast cos\theta \ast 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_{avg}$

fms中E(μ)项已经预计算出来了，还缺$E_{avg}$，也是预计算出来，不过$E_{avg}$是一维的，放到图片里存储，只能让每一行的值相同

//Eavg_MC.cpp

Vec3f IntegrateEmu(Vec3f V, float roughness, float NdotV, Vec3f Ei) {
    return Ei * NdotV * 2.0f;
}

注意前面E的图中，x轴是rough，y轴是$cos\theta$

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