环境光渲染(二）

3. 高光部分

3.1 原理

光照方程中的高光部分

L_{o s} (\vec{ω_{o}}) = \int_{Ω} f_{s} (\vec{ω_{i}}, \vec{ω_{o}}) L_{i} (\vec{ω_{i}}) ({\vec{ω}}_{i} \cdot \vec{n}) d ω_{i}

同样使用蒙特卡洛进行近似积分

\begin{aligned} \int_{Ω} f_{s} (\vec{ω_{i}}, \vec{ω_{o}}) L_{i} (\vec{ω_{i}}) ({\vec{ω}}_{i} \vec{n}) d ω_{i} \\ \approx & \frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) L_{i} (\vec{ω_{k}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})} \end{aligned}

EPIC在此基础上提出了一种近似计算方法

\begin{aligned} \frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) L_{i} (\vec{ω_{k}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}}, \vec{ω_{o}})} \\ \approx & (\frac{1}{N} \sum_{k = 1}^{N} L_{i} (\vec{ω_{k}})) (\frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})}) \end{aligned}

这种方法的好处是可以通过分开预计算，尽量减少实际渲染时的计算量，这个公式可以理解为“Env. lighting x BRDF”，左边的第一部分是小平面受到的所有环境光的采样值，第二部分只跟材质的BRDF属性相关。

3.2 预渲染环境光贴图（pre-filtered environment map)

第一个拆分项只和环境光照以及采样点有关，首先采样点集中在 ${\vec{ω}}_{r}$ 周围（ ${\vec{ω}}_{r}$ 是 ${\vec{ω}}_{o}$ 相对于法线 $\vec{n}$ 的反射向量），并且粗糙度越大，采样点越分散。实际应用中，这个拆分项是通过把环境光贴图预先渲染到一个环境贴图中实现的，称为为"pre-filter map"或者"radiance map"，最终渲染时根据 ${\vec{ω}}_{r}$ ，以及把粗糙度(0~1)转换为mip获得该方向的数值

\frac{1}{N} \sum_{k = 1}^{N} L_{i} (\vec{ω_{k}}) = TexCubeLod (PrefilterMap, {\vec{ω}}_{r}, mip)

3.3 预渲染贴图的生成

生成这张贴图时，首先假定 $\vec{n} = \vec{ω_{o}} = \vec{ω_{r}}$ ，也就是摄像机从法线方向垂直观察小平面，然后开始采样

再最终渲染时，当 $\vec{n}$ 和 $\vec{ω_{o}}$ 不等时，采样点的分布会有一定差异，但基本上可以忽略

3.3.1 采样点的生成

根据前面章节的分析，微表面的法线分布可以用法线分布函数 $D (h)$ 描述，并且符合

\int_{Ω} D (h) \cos θ_{h} d ω = 1

由于微表面法线的概率密度函数 $pdf (h)$ 符合以下公式

\int_{Ω} pdf (h) = 1

所以

pdf (h) = D (h) \cos θ_{h}

用极坐标表示

pdf (θ, ϕ) = D (h) \cos θ \sin θ

对于各向同性系统， $pdf (θ, ϕ)$ 和 $ϕ$ 无关，所以

\begin{aligned} pdf (ϕ) & = \frac{1}{2 π} \\ pdf (θ) & = \int_{0}^{2 π} pdf (θ, ϕ) d ϕ = 2 π \cdot pdf (θ, ϕ) \end{aligned}

D函数使用GGX函数，也就是

D (h) = \frac{α^{2}}{π (\cos^{2} θ (α^{2} - 1) + 1)^{2}}

其中 $α = {roughness}^{2}$ ，此时

pdf (θ) = 2 π \cdot pdf (θ, ϕ) = \frac{2 α^{2} \cos θ \sin θ}{(\cos^{2} θ (α^{2} - 1) + 1)^{2}}

计算 $θ$ 的概率分布函数

\begin{aligned} cdf (θ_{h}) & = \int_{0}^{θ_{h}} pdf (θ) d θ \\ = \int_{0}^{θ_{h}} \frac{2 α^{2} \cos θ \sin θ}{(\cos^{2} θ (α^{2} - 1) + 1)^{2}} d θ \end{aligned}

使用Mathmatica计算可得

\begin{aligned} cdf (θ) & = \frac{\sin^{2} θ}{\cos^{2} θ (α^{2} - 1) + 1} \\ = \frac{1 - \cos^{2} θ}{\cos^{2} θ (α^{2} - 1) + 1} \end{aligned}

根据前面的概率知识，可以得到微表面的法线的球坐标 $θ$ 和 $ϕ$ 的采样函数为

\begin{aligned} θ & = {cdf}^{- 1} (ξ_{θ}) = \arccos (\sqrt{\frac{1 - ξ_{θ}}{(α^{2} - 1) ξ_{θ} + 1}}) \\ ϕ & = 2 π ξ_{ϕ} \end{aligned}

然后即可计算出光线采样点的入射向量 $\vec{ω_{k}}$

3.3.2 代码

实际运算时，每个采样点的权重是不一样的，采用和宏法线的点积作为权重

\frac{1}{N} \sum_{k = 1}^{N} L_{i} ({\vec{ω}}_{k}) \approx \frac{1}{\sum_{k = 1}^{N} (\vec{n} \cdot {\vec{ω}}_{k})} \sum_{k = 1}^{N} L_{i} ({\vec{ω}}_{k}) (\vec{n} \cdot {\vec{ω}}_{k})

// make the simplyfying assumption that V equals R equals the normal 
vec3 R = N;
vec3 V = R;

const uint SAMPLE_COUNT = 1024u;
vec3 prefilteredColor = vec3_splat(0.0);
float totalWeight = 0.0;

for(uint i = 0u; i < SAMPLE_COUNT; ++i)
{
    // generates a sample vector that's biased towards the preferred 
    //alignment direction (importance sampling).
    vec2 Xi = Hammersley(i, SAMPLE_COUNT);
    vec3 H = ImportanceSampleGGX(Xi, N, roughness);
    vec3 L  = normalize(2.0 * dot(V, H) * H - V);

    float NdotL = max(dot(N, L), 0.0);
    if(NdotL > 0.0)
    {
        prefilteredColor += texture2D(s_texSkybox, SampleSphericalMap(L)).rgb * NdotL;
        totalWeight      += NdotL;
    }
}

prefilteredColor = prefilteredColor / totalWeight;

3.3 BRDF部分

观察公式的第二个拆分项

\frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})}

其中

\begin{aligned} f_{s} ({\vec{ω}}_{k}, {\vec{ω}}_{o}) & = \frac{F ({\vec{ω}}_{k}, {\vec{ω}}_{h}) D ({\vec{ω}}_{h}) G ({\vec{ω}}_{k}, {\vec{ω}}_{o})}{4 (\vec{ω_{k}} \cdot \vec{n}) (\vec{ω_{o}} \cdot \vec{n})} \\ pdf (\vec{ω_{k}}) & = \frac{pdf (\vec{ω_{h}})}{4 (\vec{ω_{h}} \cdot \vec{ω_{o}})} = \frac{D (\vec{ω_{h}}) (\vec{ω_{h}} \cdot \vec{n})}{4 (\vec{ω_{h}} \cdot \vec{ω_{o}})} \\ F ({\vec{ω}}_{k}, {\vec{ω}}_{h}) & = F_{0} + (1 - F_{0}) (1 - \vec{ω_{h}} \cdot \vec{ω_{o}})^{5} \end{aligned}

设 $F_{c} = (1 - \vec{ω_{h}} \cdot \vec{ω_{o}})^{5}$ ，那么

F ({\vec{ω}}_{k}, {\vec{ω}}_{h}) = F_{0} (1 - F_{c}) + F_{c}

带入上面的式中，可得

\begin{aligned} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})} & = F_{0} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}})}{F ({\vec{ω}}_{k}, {\vec{ω}}_{h}) pdf (\vec{ω_{k}})} (1 - F_{c}) (\vec{ω_{k}} \cdot \vec{n}) + \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}})}{F ({\vec{ω}}_{k}, {\vec{ω}}_{h}) pdf (\vec{ω_{k}})} F_{c} (\vec{ω_{k}} \cdot \vec{n}) \\ = F_{0} (1 - F_{c}) \frac{G ({\vec{ω}}_{k}, {\vec{ω}}_{o}) (\vec{ω_{h}} \cdot \vec{ω_{o}})}{(\vec{ω_{h}} \cdot \vec{n}) (\vec{ω_{o}} \cdot \vec{n})} + F_{c} \frac{G ({\vec{ω}}_{k}, {\vec{ω}}_{o}) (\vec{ω_{h}} \cdot \vec{ω_{o}})}{(\vec{ω_{h}} \cdot \vec{n}) (\vec{ω_{o}} \cdot \vec{n})} \end{aligned}

设

G_{v} = \frac{G ({\vec{ω}}_{k}, {\vec{ω}}_{o}) (\vec{ω_{h}} \cdot \vec{ω_{o}})}{(\vec{ω_{h}} \cdot \vec{n}) (\vec{ω_{o}} \cdot \vec{n})}

那么可得

\frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})} = F_{0} (\frac{1}{N} \sum (1 - F_{c}) G_{v}) + \frac{1}{N} \sum F_{c} G_{v}

3.4 BRDF查找贴图

上面的公式中，除了 $F_{0}$ 和运行时的材质相关，其他数值都只取决于所使用的微表面的函数类型，所以可以预先计算在一张二维贴图中

在运行时可以使用查表得方式从这张贴图上获取数据

\frac{1}{N} \sum_{k = 1}^{N} \frac{f_{s} (\vec{ω_{k}}, \vec{ω_{o}}) (\vec{ω_{k}} \cdot \vec{n})}{pdf (\vec{ω_{k}})} = F_{0} * Tex2D (\vec{ω_{o}} \cdot \vec{n}, Roughness) . r + Tex2D (\vec{ω_{o}} \cdot \vec{n}, Roughness) . g

3.4.1 BRDF贴图的生成

vec2 IntegrateBRDF(float NdotV, float roughness)
{
    vec3 V;
    V.x = sqrt(1.0 - NdotV*NdotV);
    V.y = 0.0;
    V.z = NdotV;

    float A = 0.0;
    float B = 0.0; 

    vec3 N = vec3(0.0, 0.0, 1.0);
    
    const uint SAMPLE_COUNT = 1024u;
    for(uint i = 0u; i < SAMPLE_COUNT; ++i)
    {
        // generates a sample vector that's biased towards the
        // preferred alignment direction (importance sampling).
        vec2 Xi = Hammersley(i, SAMPLE_COUNT);
        vec3 H = ImportanceSampleGGX(Xi, N, roughness);
        vec3 L = normalize(2.0 * dot(V, H) * H - V);

        float NdotL = max(L.z, 0.0);
        float NdotH = max(H.z, 0.0);
        float VdotH = max(dot(V, H), 0.0);

        if(NdotL > 0.0)
        {
            float G = GeometrySmith(N, V, L, roughness);
            float G_Vis = (G * VdotH) / (NdotH * NdotV);
            float Fc = pow(1.0 - VdotH, 5.0);

            A += (1.0 - Fc) * G_Vis;
            B += Fc * G_Vis;
        }
    }
    A /= float(SAMPLE_COUNT);
    B /= float(SAMPLE_COUNT);
    return vec2(A, B);
}