古法渲染：DisneyBSDF的实现

这是一篇 DisneyBSDF 实现笔记，主要记录踩坑过程中的三个核心模块：

• 玻璃材质BSDF的实现细节；

• 玻璃材质BSDF 与 DisneyBRDF 的整合；

• VNDF 可见性法线分布采样。

在玻璃材质 BSDF 分支中，个人感觉最棘手的部分是 透射率（transmittance） 的处理。由于光线在介质边界来回折射，参数和方向需要反复转换，实现时很容易搞混符号和相对关系，写起来颇为绕手。

将玻璃材质 BSDF 整合到 DisneyBRDF 中则相对直接：按一定概率将玻璃材质 BSDF 采样策略与原有 BRDF 采样策略做混合（例如根据表面透射率进行 MIS 或简单分支选择）即可。

VNDF（Visible Normal Distribution Function）采样是一种降低蒙特卡洛方差的重要技巧。其实现本身并不复杂，但由于采样分布从普通 NDF 变为可见性加权分布，对应的 CDF/PDF 也发生了变化。因此在蒙特卡洛积分时，必须根据新的采样分布重新计算 PDF，才能保证估计量的无偏性。

不过话又说回来，虽然我采用了 VNDF 采样，但实际渲染结果好像并没有多大变化——至少我自己没看出来。不确定是哪里实现得有问题，还是效果本来就是这样。目前也不打算深究了，以后遇到问题再说。

另外，本文不会涉及 DisneyBRDF 本身的实现 。Disney 官方很久之前就已经开源了源码，你可以在这里找到这份经典的代码：DisneyBRDF

网络上也有不少优秀的实现案例，本文很大程度上借鉴了 EzRT 的版本：EzRT/DisneyBRDF。

OK，闲话到此位置，先放一张渲染图~

玻璃材质 BSDF 分支

玻璃材质 BSDF 分支 是实现 Disney BSDF 的关键部分。不过透射过程通常并不是独立发生的。对于玻璃这类透明材质，光线在与表面发生作用时，往往会同时出现反射与透射两种行为。因此，在讨论 BTDF 的实现时，我打算一并介绍玻璃材质下 BTDF 与 BRDF （即 BSDF）的计算与采样过程。

这一小节主要分为以下三个部分：

1. 玻璃材质下 BTDF 与 BRDF 评估
2. 玻璃材质下 BTDF 与 BRDF 的采样方法
3. 玻璃材质下 BTDF 与 BRDF 的概率密度计算

为了更准确地描述各个参数之间的关系，先约定本文使用的符号：

• $n_i$：入射介质的折射率；
• $n_o$：出射介质的折射率；
• $\eta$: 入射介质与出射介质折射率之比，即 $\frac{n_i}{n_o}$;
• $\mathbf v$: 归一化的三维向量；
• $\mathbf m$: 微平面法线方向；
• $\mathbf h$: 采样到的半程向量；
• $\mathbf n$: 宏观平面法线方向；
• $\mathbf i$：入射方向，从表面交点指向光源，也可记作 $L$ 或 $\boldsymbol \omega_i$；
• $\mathbf o$：出射方向，从表面交点指向视点，也可记作 $V$ 或 $\boldsymbol \omega_o$；

这里我为入射方向和出射方向约定的不同的符号，这是因为后面两个符号用在代码中会更加清楚一些，而前面的$\mathbf i$ 和 $\mathbf o$ 则更加适合数学表达（好吧我承认有些多余）。

一个简单的示意图如下：

示意图中，上方介质的折射率为$n_1$，下方介质的折射率为 $n_2$。

玻璃材质下 BTDF 与 BRDF 评估

在之前的文章中，我们已经详细讨论了 BTDF 与 BRDF 的具体形式，这里只做简单回顾。

基于微表面模型评估 BSDF（包括 BRDF 和 BTDF）通常需要以下三个组成部分：

1. 法线分布函数（NDF） $D$ ；

2. 几何遮蔽（阴影-遮蔽）函数 $G$ ；

3. 菲涅尔项 $F$ 。

法线分布函数 $D$

法线分布函数采用 GGX（Trowbridge-Reitz）分布，即 Disney BRDF 论文附录中 GTR 分布取 γ=2 的情形。这里先不考虑各向异性，其公式为：

$$D_{\mathrm{GTR}_2}(\mathbf n, \mathbf h,\alpha)=\frac{\alpha^2}{\pi(1+(\alpha^2-1)|\mathbf n \cdot \mathbf h|^2)^2}$$

其中，$\mathbf h$ 为半程向量（Half-vector），$\mathbf n$ 为宏观表面法线，$\alpha$ 为与粗糙度相关的参数，通常取

$$\alpha = roughness^2$$

对应的代码实现如下：：

float GTR2(float NdotH, float a)
{
    a = max(a, 1e-3);
    float a2 = a*a;
    float t = 1 + (a2-1)*NdotH*NdotH;
    return a2 / (PI * t*t);
}

阴影-遮蔽函数 $G$

阴影-遮蔽函数采用 Smith GGX 模型。其中，$G_1$​ 为单方向遮蔽函数（Masking Function），表示从方向 $\mathbf v$ 观察时，微表面未被几何遮蔽的可见比例：

$$G_1(\mathbf v,\,\alpha)=\frac{2|\mathbf v \cdot \mathbf n|}{|\mathbf v \cdot \mathbf n|+\sqrt{\alpha^2+(1-\alpha^2)|\mathbf v \cdot \mathbf n|^2}}$$

而 $G_2$​ 为联合阴影-遮蔽函数，描述光线在入射和出射过程中均未被遮挡的概率。在 Smith 近似下，入射侧与出射侧的遮蔽事件被视为相互独立，因此 $G_2$​ 可分解为两个 $G_1$​ 的乘积：

$$G_2(\mathbf i, \,\mathbf o, \,\alpha)=G_1(\mathbf i,\,\alpha)\cdot G_1(\mathbf o, \,\alpha)$$

对应的代码实现如下：

float smithG1_GGX(float NdotV, float alpha)
{
    NdotV = max(NdotV, 1e-6);
    alpha = max(alpha, 1e-3);

    float a2 = alpha * alpha;
    float n2 = NdotV * NdotV;
    return (2.0 * NdotV) / (NdotV + sqrt(a2 + (1.0 - a2) * n2));
}

float smithG_GGX_Glass(float NdotL, float NdotV, float alpha)
{
    return smithG1_GGX(NdotL, alpha) * smithG1_GGX(NdotV, alpha);
}

在 DisneyBRDF 源码中，为了简化计算，阴影-遮蔽项与 BRDF 分母中的余弦项被合并打包成了一个辅助函数（即 visibility 项），使得评估时可以直接写成 $D \cdot vis \cdot F$ 的紧凑形式。不过为了追求逻辑清晰并确保正确性，我还是选择显式实现完整的 $G_2$​ ，并在 BSDF 评估中使用标准公式。（这玩意坑了我一整天才找到这个bug）

菲涅尔项 F

在实时渲染中，菲涅尔项通常使用 Schlick 近似来简化计算。然而，当光线从光密介质进入光疏介质（即从高折射率透射到低折射率）时，若入射角大于临界角，会发生全反射（TIR）。此时近似公式不再准确，甚至可能给出错误的透射率，导致能量不守恒。

虽然可以通过额外的判断和修正来处理这种情况，但为了确保反射与透射分支中能量分配的正确性，我们直接使用完整的 Fresnel 公式进行计算。其表达式为：

$$\begin{align*} & F_{g}=\frac{1}{2}\left(R_s^2+R_p^2\right) \\ & R_{s}=\frac{\mathbf i \cdot \mathbf h -\eta \,(\mathbf o \cdot \mathbf h)}{\mathbf i \cdot \mathbf h+\eta\,(\mathbf o \cdot \mathbf h)} \\ & R_{p}=\frac{\eta \,(\mathbf i \cdot \mathbf h)- \mathbf o \cdot \mathbf h}{\eta \,(\mathbf i \cdot \mathbf h)+\mathbf o \cdot \mathbf h} \end{align*}$$

对应的代码实现如下：：

float FrDielectric(float cosThetaI, float eta)
{
    cosThetaI = clamp(cosThetaI, -1.0, 1.0);

    float sin2ThetaI = max(0.0, 1.0 - cosThetaI * cosThetaI);
    float sin2ThetaT = eta * eta * sin2ThetaI;

    if (sin2ThetaT >= 1.0)
        return 1.0;

    float cosThetaT = sqrt(max(0.0, 1.0 - sin2ThetaT));
    float absCosI = abs(cosThetaI);

    float Rs = (absCosI - eta * cosThetaT) / max(absCosI + eta * cosThetaT, 1e-6);
    float Rp = (eta * absCosI - cosThetaT) / max(eta * absCosI + cosThetaT, 1e-6);

    return 0.5 * (Rs * Rs + Rp * Rp);
}

玻璃材质的 BRDF 评估函数

对于反射路径，Disney 玻璃材质的 BRDF 采用标准微表面镜面反射模型，评估公式如下：

$$f_r(\mathbf{i},\mathbf{o},\mathbf{n})=\mathrm{BaseColor}\left(\frac{F(\mathbf{i},\mathbf{h_r})G(\mathbf{i},\mathbf{o},\mathbf{h_r})D(\mathbf{h_r})}{4\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{o}\cdot\mathbf{n}\right|}\right)$$

其中 $\mathbf{h_{r}}=\frac{\mathbf{i}+\mathbf{o}}{|\mathbf{i}+\mathbf{o}|}$​ 为反射半程向量，$F$ 为菲涅尔反射系数，$G$ 为联合阴影-遮蔽函数，$D$ 为法线分布函数，α 为粗糙度参数。分母中的 $4|\mathbf{i}\cdot\mathbf{n}||\mathbf{o}\cdot\mathbf{n}|$ 是微表面 BRDF 的标准归一化项。

对应的代码实现如下：

vec3 HalfVectorReflection(vec3 wi, vec3 wo)
{
    vec3 h = wi + wo;
    float len2 = dot(h, h);
    if (len2 <= 1e-12)
        return vec3(0.0);
    return h * inversesqrt(len2);
}

vec3 EvalGlassReflection(vec3 wi, vec3 wo, vec3 Ns, float eta, int MaterialSlotID)
{
    if (!SameHemisphere(wi, wo, Ns))
        return vec3(0.0);

    float NdotL = AbsDot(Ns, wi);
    float NdotV = AbsDot(Ns, wo);
    if (NdotL <= 1e-6 || NdotV <= 1e-6)
        return vec3(0.0);

    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float roughness = max(Mat.Roughness, 0.01);
    float alpha = roughness * roughness;

    vec3 H = HalfVectorReflection(wi, wo);
    if (dot(H, H) <= 0.0)
        return vec3(0.0);

    if (dot(H, Ns) < 0.0)
        H = -H;

    float NdotH = AbsDot(Ns, H);
    float VdotH = AbsDot(wo, H);
    if (NdotH <= 1e-6 || VdotH <= 1e-6)
        return vec3(0.0);

    float Fg = FrDielectric(dot(wi, H), eta);
    float Dg = GTR2(NdotH, alpha);
    float Gg = smithG_GGX_Glass(NdotL, NdotV, alpha);

    float denom = 4.0 * NdotL * NdotV;

    return Mat.BaseColor * Fg * Dg * Gg / denom;
}

玻璃材质的 BTDF 评估函数

对于透射（折射）路径，Disney 玻璃材质的 BTDF 评估公式如下。菲涅尔项使用透射部分 $(1−F)$ ，而 BaseColor 取平方根以近似模拟光线在介质内部往返时的体积吸收（Beer 定律的简化处理）：

$$f_{t}\left(\mathbf{i},\mathbf{o},\mathbf{n}\right)=\sqrt{\mathrm{BaseColor}}\frac{\left|\mathbf{i}\cdot\mathbf{h}_{\mathrm{t}}\right|\left|\mathbf{o}\cdot\mathbf{h}_{\mathrm{t}}\right|}{\left|\mathbf{i}\cdot\mathbf{n}\right|\left|\mathbf{o}\cdot\mathbf{n}\right|}\frac{\left(1-F(\mathbf{i},\mathbf{h}_{\mathrm{t}})\right)G(\mathbf{i},\mathbf{o},\mathbf{h}_{\mathrm{t}})D(\mathbf{h}_{\mathrm{t}})}{\left(\eta\,(\mathbf{i}\cdot\mathbf{h}_{\mathrm{t}})+(\mathbf{o}\cdot\mathbf{h}_{\mathrm{t}})\right)^{2}}$$

其中 $\mathbf{h_{t}}$​ 为折射半程向量（即满足 Snell 定律的微表面法线方向），$η=\frac{n_i}{n_o}$​ 为相对折射率，$G$ 与 $D$ 的形式与反射分支相同，均由粗糙度 $\alpha$ 控制。

写成代码：

bool SameHemisphere(vec3 a, vec3 b, vec3 N)
{
    return dot(a, N) * dot(b, N) > 0.0;
}

vec3 HalfVectorTransmission(vec3 wi, vec3 wo, float eta_wi2wo)
{
    vec3 h =  wi * eta_wi2wo + wo;
    float len2 = dot(h, h);
    if (len2 <= 1e-12)
        return vec3(0.0);
    return h * inversesqrt(len2);
}

vec3 EvalGlassTransmission(vec3 wi, vec3 wo, vec3 Ns, float eta_wi2wo, int MaterialSlotID)
{
    if (SameHemisphere(wi, wo, Ns))
        return vec3(0.0);

    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float NdotL = AbsDot(Ns, wi);
    float NdotV = AbsDot(Ns, wo);
    if (NdotL <= 1e-6 || NdotV <= 1e-6)
        return vec3(0.0);

    float roughness = max(Mat.Roughness, 0.01);
    float alpha = roughness * roughness;

    vec3 H = HalfVectorTransmission(wi, wo, eta_wi2wo);
    if (dot(H, H) <= 0.0)
        return vec3(0.0);

    if (dot(H, Ns) < 0.0)
        H = -H;

    float cosWi = dot(Ns, wi);
    float cosWo = dot(Ns, wo);
    if (dot(H, wi) * cosWi < 0.0 || dot(H, wo) * cosWo < 0.0)
        return vec3(0.0);

    float NdotH = AbsDot(Ns, H);
    float VdotH = AbsDot(wo, H);
    float LdotH = AbsDot(wi, H);

    if (NdotH <= 1e-6 || VdotH <= 1e-6 || LdotH <= 1e-6)
        return vec3(0.0);

    float Fg = FrDielectric(dot(wi, H), eta_wi2wo);
    float Dg = GTR2(NdotH, alpha);
    float Gg = smithG_GGX_Glass(NdotL, NdotV, alpha);


    float sqrtDenom = eta_wi2wo * LdotH + VdotH;
    float denom = max(NdotL * NdotV * sqrtDenom * sqrtDenom, 1e-6);

    vec3 tint = sqrt(max(Mat.BaseColor, vec3(0.0)));

    return tint * abs((1.0 - Fg) * Dg * Gg * LdotH * VdotH / denom);
}

玻璃材质的 BSDF 评估函数

现在要做的就是将玻璃材质的 BTDF 和 BRDF 整合成 BSDF 评估函数.

在路径追踪中，一次 bounce 只会对应一条被采样出来的出射方向，因此在实际评估时，我们需要先判断当前的 $\mathbf i$ 和 $\mathbf o$ 究竟对应的是反射事件还是透射事件，然后再进入对应的分支进行计算。这个判断最直接的方法，就是看入射方向和出射方向是否位于几何法线 $N_g$​ 的同一侧。代码如下：

vec3 BSDF(vec3 wi, vec3 wo, vec3 Ns, vec3 Ng, bool bIsInside, int MaterialSlotID)
{
    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float eta_wi2wo = bIsInside ?  (1.0 / max(Mat.IOR, 1.0001)) : max(Mat.IOR, 1.0001);
    bool sameSide = dot(Ng, wi) * dot(Ng, wo) > 0.0;

    vec3 f = vec3(0.0);

    if (sameSide)
        f += EvalGlassReflection(wi, wo, Ns, 1.0 / eta_wi2wo, MaterialSlotID);
    else
        f += EvalGlassTransmission(wi, wo, Ns, eta_wi2wo, MaterialSlotID);

    return f;
}

这里有一个非常容易出错、但又非常关键的细节：反射分支和透射分支传入的 $\eta$ 并不相同。

原因在于，玻璃材质的 BSDF 在计算 Fresnel 项时，所使用的折射率比需要根据当前对应的是反射光线还是折射光线来决定。虽然看起来只是“传参顺序不同”，但如果这里处理错了，最终结果往往会出现明显偏差，而且调试起来也很麻烦。

还是这张示意图：

在这个示意关系下：

• 如果 $\mathbf{i}$ 对应的是反射方向，那么 Fresnel 计算中使用的折射率比应为

  $$\eta = \frac{n_1}{n_2}$$

• 如果 $\mathbf{i}$ 对应的是折射方向，那么应使用

  $$\eta = \frac{n_2}{n_1}$$

换句话说，同样是一组入射 / 出射方向，在进入不同分支时，$\eta$ 的定义方向也要随之切换。这一点在实现玻璃 BSDF 时尤其重要，也是很容易出错的地方。

玻璃材质下 BTDF 与 BRDF 的采样方法

在玻璃材质中，BTDF（透射） 和 BRDF（反射） 的采样仍然沿用比较经典的思路：先对半程向量 $\mathbf h$ 进行重要性采样，再根据 $\mathbf h$ 与已知方向之间的几何关系，反推出对应的入射方向 $\mathbf i$。

整个过程可以概括为下面几个步骤：

1. 重要性采样半程向量 $\mathbf i$

2. 路径追踪是一个逆向构造光路的过程。也就是说，我们是从视点出发，逐步回溯一条可能通向光源的路径。因此在采样阶段，已知量通常是出射方向 $\mathbf o$；

3. 根据 $\mathbf h$ 与 $\mathbf o$ 的夹角关系计算 Fresnel 项。有了 Fresnel 系数之后，就可以进一步得到本次采样中选择反射还是折射的概率；

4. 生成一个随机数，并按照对应概率在反射分支和折射分支之间进行选择。

这里有一个需要特别注意的细节：我们现在是已知出射方向 $\mathbf o$，反推出入射方向 $\mathbf i$。因此，在折射分支中使用的折射率比，不能再直接照搬“入射到出射”的写法，而要改成从 $\mathbf o$ 所在介质到 $\mathbf i$ 所在介质的折射率之比，即

$$\eta_{\boldsymbol\omega_o \rightarrow \boldsymbol\omega_i} = \frac{n_o}{n_i}$$

这个方向性问题并不仅仅只是只是符号上的差别。一旦这里的折射率比写反，后续的折射方向计算、Fresnel 概率，甚至整个采样分支的结果都会出错。

相关代码如下：

struct SampleBSDFResult
{
    vec3 wi;
    float cosTheta;
    float PDF;
};

SampleBSDFResult SampleBSDF_Glass(vec3 wo, vec3 Ns, bool bIsInside, int MaterialSlotID)
{
    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float roughness = max(Mat.Roughness, 0.01);
    float alpha = roughness * roughness;
    float eta_wo2wi = bIsInside ? max(Mat.IOR, 1.0001) : (1.0 / max(Mat.IOR, 1.0001));

    SampleBSDFResult result;
    result.wi = vec3(0.0);
    result.cosTheta = 0.0;
    result.PDF = 0.0;

    vec3 H = SampleGTR2(rand(), rand(), Ns, alpha);

    if (dot(wo, H) < 0.0)
        H = -H;

    float Fg = FrDielectric(dot(wo, H), eta_wo2wi);

    vec3 wi_reflect = reflect(-wo, H);
    vec3 wi_refract = refract(-wo, H, eta_wo2wi);

    bool tir = dot(wi_refract, wi_refract) <= 1e-12; // 全反射
    float xi = rand();

    if (tir || xi < Fg)
    {
        result.wi = wi_reflect;

        if (dot(Ns, result.wi) <= 0.0)
        {
            result.wi = vec3(0.0);
            return result;
        }
    }
    else
    {
        result.wi = wi_refract;

        if (dot(Ns, result.wi) * dot(Ns, wo) >= 0.0)
        {
            result.wi = vec3(0.0);
            return result;
        }
    }

    result.cosTheta = AbsDot(Ns, result.wi);
    return result;
}

玻璃材质下 BTDF 与 BRDF 的概率密度计算

在采样阶段，我们首先是根据法线分布函数 $D$ 对半程向量 $\mathbf h$ 进行采样的。借助微表面模型，这个分布函数在投影到表面面积后天然满足归一化条件

$$\int_\Omega D(\mathbf m) |\mathbf n \cdot \mathbf m| \,\mathrm d \omega = 1$$

因此，采样到某个半程向量 $\mathbf h$ 的概率密度可以写成：

$$p(\mathbf h) = D(\mathbf h)\, |\mathbf n \cdot \mathbf h|$$

不过，这还不是路径追踪真正需要的量。对于 BSDF 采样来说，我们最终关心的是入射方向 $\mathbf i$ 的概率密度，而不是半程向量 $\mathbf h$ 的概率密度。因此，还需要借助一个雅可比项，把从 $\mathbf h$ 空间得到的概率密度转换到 $\mathbf i$ 空间：

$$p(\mathbf i) = p(\mathbf h)\left\|\frac{\partial \boldsymbol \omega_h}{\partial \boldsymbol \omega_i}\right\|$$

这个转换过程在反射分支和折射分支中并不相同，因此两者对应的雅可比项也不同。

对于反射分支，雅可比项为：

$$\left\|\frac{\partial\omega_{\mathbf{h_r}}}{\partial\omega_{\mathbf{i}}}\right\| = \frac{1}{4\left|\mathbf{i} \cdot\mathbf{h_r}\right|}$$

对于折射分支，雅可比项为：

$$\left\|\frac{\partial\omega_{\mathbf{h_t}}}{\partial\omega_{\mathbf{i}}}\right\| = \frac{\eta_{i}^{2}\left|\mathbf{i}\cdot\mathbf{h_{t}}\right|}{\left(\eta_{o}(\mathbf{o}\cdot\mathbf{h_{t}})+\eta_{i}(\mathbf{i}\cdot\mathbf{h_{t}})\right)^{2}}$$

除了雅可比项之外，还有一个细节也不能忽略：在前面的采样步骤里，我们是利用 Fresnel 项 来决定当前样本落到反射分支还是折射分支的。换句话说，Fresnel 系数本身也是这两个分支的离散选择概率，因此在计算最终 PDF 时，还需要把这部分概率一并乘进去。

于是，玻璃材质下反射分支的概率密度可以写成：

$$p_r(\mathbf i) = F \cdot p(\mathbf h)\left\|\frac{\partial\omega_{\mathbf{h_r}}}{\partial\omega_{\mathbf{i}}}\right\|$$

而折射分支的概率密度为：

$$p_t(\mathbf i) = (1 - F)\cdot p(\mathbf h)\left\|\frac{\partial\omega_{\mathbf{h_t}}}{\partial\omega_{\mathbf{i}}}\right\|$$

把这些环节串起来之后，整个玻璃材质的采样 PDF 才算完整闭合。

将这部分写成代码如下：

float PdfGlassReflectionIS(vec3 wi, vec3 wo, vec3 Ns, float eta_wo2wi, float alpha)
{
    if (!SameHemisphere(wi, wo, Ns))
        return 0.0;

    vec3 H = HalfVectorReflection(wi, wo);
    if (dot(H, H) <= 0.0)
        return 0.0;

    if (dot(H, Ns) < 0.0)
        H = -H;

    float NdotH = AbsDot(Ns, H);
    float VdotH = AbsDot(wo, H);

    if (NdotH <= 1e-6 || VdotH <= 1e-6)
        return 0.0;

    float pdfH = GTR2(NdotH, alpha) * NdotH;
    float Fg = FrDielectric(dot(wo, H), eta_wo2wi);

    return  Fg * pdfH / max(4.0 * VdotH, 1e-6);
}

float PdfGlassTransmissionIS(vec3 wi, vec3 wo, vec3 Ns, float eta_wo2wi, float alpha)
{
    if (SameHemisphere(wi, wo, Ns))
        return 0.0;

    float eta_wi2wo = 1.0/eta_wo2wi;
    vec3 H = HalfVectorTransmission(wi, wo, eta_wi2wo);
    if (dot(H, H) <= 0.0)
        return 0.0;

    if (dot(H, Ns) < 0.0)
        H = -H;

    float NdotH = AbsDot(Ns, H);
    float VdotH = AbsDot(wo, H);
    float LdotH = AbsDot(wi, H);

    if (NdotH <= 1e-6 || VdotH <= 1e-6 || LdotH <= 1e-6)
        return 0.0;

    float pdfH = GTR2(NdotH, alpha) * NdotH;
    float Fg = FrDielectric(dot(wo, H), eta_wo2wi);

    float sqrtDenom = eta_wo2wi * VdotH + LdotH;
    float dwh_dwi = abs((LdotH) / max(sqrtDenom * sqrtDenom, 1e-6));

    return (1 - Fg) * pdfH * dwh_dwi;
}

float BSDF_Glass_PDF_Eval(vec3 wi, vec3 wo, vec3 Ns, bool bIsInside, int MaterialSlotID)
{
    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float roughness = max(Mat.Roughness, 0.01);
    float alpha = roughness * roughness;
    float eta_wo2wi = bIsInside ? max(Mat.IOR, 1.0001) : (1.0 / max(Mat.IOR, 1.0001));

    if (SameHemisphere(wi, wo, Ns))
        return PdfGlassReflectionIS(wi, wo, Ns, eta_wo2wi, alpha);
    else
        return PdfGlassTransmissionIS(wi, wo, Ns, eta_wo2wi, alpha);
}

这种采样方式被称为 One-Sample MIS，也就是单样本多重重要性采样。它的核心思想是：虽然我们手头有多种不同的采样策略，但在每次采样时只真正选择其中一种来生成样本；只要估计器和对应的权重设计正确，最终仍然可以保证估计结果是无偏的。

由于下一节还会继续用到这一部分，这里先给出一个简洁的无偏性说明。

在 One-Sample MIS 的框架下，假设共有 $K$ 种采样策略，则估计器可以写成：

$$\hat I(i, X) = \frac{w_i(X) f(X)}{c_i p_i(X)}$$

其中：

• $c_i$​ 表示选择第 $i$ 个采样策略的概率；

• $p_i(X)$ 表示策略 $i$ 对样本 $X$ 的概率密度；

• $w_i(X)$ 是对应的 MIS 权重，并满足

$$\sum_{i=1}^{K} w_i(X) = 1$$

回到前面玻璃材质的反射 / 折射例子中，这个过程其实很好理解：当采样结果落在反射方向上时，反射分支的权重为 1，折射分支的权重为 0；反过来，当采样结果落在折射方向上时，反射分支的权重为 0，折射分支的权重为 1。至于每个分支被选中的概率 $c_i$​，则正是由 Fresnel 项 决定的。

下面对这个估计器求期望：

$$\begin{align*} E[\hat I] &= \sum_{i=1}^{K} c_i \int \frac{w_i(X) f(X)}{c_i p_i(X)} p_i(X)\,\mathrm dX \\ &= \sum_{i=1}^{K} \int w_i(X) f(X)\,\mathrm dX \\ &= \int \sum_{i=1}^{K} w_i(X) f(X)\,\mathrm dX \\ &= \int f(X)\,\mathrm dX \\ &= I \end{align*}$$

因此，只要权重满足归一性条件，One-Sample MIS 构造出来的估计器就是无偏的。

更进一步地，如果采用最常见的 balance heuristic，也就是令权重取为

$$w_i(X) = \frac{c_i p_i(X)}{\sum_j c_j p_j(X)}$$

那么估计器可以进一步化简为：

$$\hat I(i, X) = \frac{f(X)}{\sum_j c_j p_j(X)}$$

此时，整个混合采样过程对应的概率密度函数也就变成了

$$p_{\mathrm{mix}}(X) = \sum_j c_j p_j(X)$$

这个形式在实现上非常方便，因为它直接把多种采样策略合并成了一个统一的混合 PDF。下一节里我们所使用的估计器，正是这种写法。

附上两张渲染图：

玻璃材质BSDF 与 DisneyBRDF 的整合

将 BTDF 分支整合进 Disney BRDF 的过程其实并不复杂。这里同样可以沿用 One-Sample MIS 的思路：把“玻璃材质 BSDF”和“Disney BRDF”进一步抽象为两种不同的采样策略，再按照一定概率在两者之间进行混合。

需要注意的是，这里的“两个策略”本身都不是单一分支。比如，玻璃材质 BSDF 内部还包含反射和折射两种采样策略；而 Disney BRDF 内部则可能包含漫反射、高光、清漆等多个子策略。不过从更高一层来看，它们仍然可以分别视为两个大的采样器，再用一次 One-Sample MIS 在二者之间做选择。

这里将玻璃材质 BSDF 的选择概率定义为：

$$c_{\mathrm{glass}} = (1 - \mathrm{metallic}) \cdot \mathrm{transmission}$$

这个表达式也比较直观：只有在材质不是金属、并且具有一定透射属性时，玻璃分支才应该参与采样。得到这个概率之后，只需要生成一个随机数 $\xi$：当 $\xi < c_{\mathrm{glass}}$​ 时，进入上一节介绍的玻璃材质 BSDF 分支；否则就进入 Disney BRDF 分支。

整合了玻璃材质 BSDF 之后，整个 Disney BRDF 的采样流程可以用下图表示：

相关代码如下：

SampleBSDFResult SampleBSDF(vec3 wo, vec3 Ns, vec3 Ng, bool bIsInside, int MaterialSlotID)
{
    SampleBSDFResult result;
    result.wi = vec3(0.0);
    result.cosTheta = 0.0;
    result.PDF = 0.0;

    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float r_diffuse   = (1.0 - Mat.Metallic) * (1.0 - Mat.Transmission);
    float r_specular  = 1.0 - Mat.Transmission * (1.0 - Mat.Metallic);
    float r_glass = (1.0 - Mat.Metallic) * Mat.Transmission;
    float r_clearcoat = 0.25 * Mat.Clearcoat;
    float r_sum       = max(r_diffuse + r_specular + r_glass + r_clearcoat, 1e-8);

    float p_glass = r_glass / r_sum;

    float xi = rand();

    if (xi < p_glass)
    {
        result = SampleBSDF_Glass(wo, Ns, bIsInside, MaterialSlotID);
    }
    else
    {
        result = SampleBRDF(wo, Ns, MaterialSlotID);
    }

    if (dot(result.wi, result.wi) <= 0.0)
        return result;

    result.cosTheta = abs(dot(Ns, result.wi));
    result.PDF = max(BSDF_PDF_Eval(result.wi, wo, Ns, p_glass, bIsInside, MaterialSlotID), 1e-8);

    return result;
}

除了采样策略需要混合之外，BSDF 的评估值本身也要以相同的方式进行组合。既然现在我们同时拥有 Glass BSDF 和 Disney BRDF 两部分贡献，那么最终的材质响应也可以按照相同的混合因子进行线性插值（这里是一个可替换的策略）。这里仍然直接使用 $c_{\mathrm{glass}}$​ 作为混合系数：

$$f(X) = c_{\mathrm{glass}}\cdot f_{\mathrm{glass}}(X)+(1-c_{\mathrm{glass}})f_r(X)$$

这样一来，采样策略和函数评估就在同一套概率模型下保持了一致：一方面用 $c_{\mathrm{glass}}$​ 控制样本来自玻璃分支还是 Disney BRDF 分支，另一方面也用同样的系数对两部分的 BSDF 值进行混合。这样处理后，也很容易与前面 One-Sample MIS 的推导对应起来。

相关代码如下：

vec3 BSDF(vec3 wi, vec3 wo, vec3 Ns, vec3 Ng, float p_glass, bool bIsInside, int MaterialSlotID)
{
    BSDFMaterial Mat = DecodeBSDFMaterial(MaterialSlotID);

    float eta_wi2wo = bIsInside ?  (1.0 / max(Mat.IOR, 1.0001)) : max(Mat.IOR, 1.0001);
    bool sameSide = dot(Ng, wi) * dot(Ng, wo) > 0.0;

    vec3 f = vec3(0.0);

    f += (1.0 - p_glass) * EvalBRDF(wi, wo, Ns, MaterialSlotID);

    if (sameSide)
        f += p_glass * EvalGlassReflection(wi, wo, Ns, 1.0 / eta_wi2wo, MaterialSlotID);
    else
        f += p_glass * EvalGlassTransmission(wi, wo, Ns, eta_wi2wo, MaterialSlotID);

    return f;
}

VNDF 可见性法线分布采样

前面几节里，我们一直是在用 NDF（法线分布函数） 对半程向量进行采样。这样的做法在大多数情况下没有问题，但它有一个比较明显的缺点：当视线以很低的掠射角看向表面时，会有相当一部分微表面实际上已经被遮蔽，根本不可能对当前观察方向产生贡献。

这意味着，如果在这种情况下仍然直接按照 NDF 采样，就会频繁采到那些实际上不可见的微表面法线，从而引入较大的采样方差。为了缓解这个问题，采样时就不能只考虑 NDF 本身，还需要把单向阴影遮蔽函数 $G_1$ 一并纳入进来。

如图所示，黑色区域表示那些对当前视线不可见的微平面法线。既然这些法线不会产生有效贡献，那么在采样时也就应当尽量避开它们。

这正是 VNDF（Visible Normal Distribution Function，可见性法线分布函数）采样 的核心思想：不再对所有可能的微表面法线一视同仁，而是直接对“当前视线下可见的那部分法线分布”进行采样。

在 VNDF 采样框架下，用于采样半程向量的 PDF 会变成：

$$D_{\mathbf{h}}(\mathbf{h})=\frac{G_1(\mathbf{o}, \alpha)\left|\mathbf o\cdot\mathbf{h}\right|D(\mathbf{h})}{\left|\mathbf{o}\cdot\mathbf n\right|},$$

这个式子的直观含义并不难理解：它在原本的法线分布 $D(\mathbf h)$ 基础上，额外引入了可见性项 $G_1$​，从而把那些已经被遮蔽的微表面法线“过滤”掉。分母中的 $\left|\mathbf{o}\cdot\mathbf n\right|$ 则用于完成归一化，使整个分布仍然是合法的概率密度函数。

论文Sampling the GGX Distribution of Visible Normals已经给出了核心采样代码，直接copy：

// Input Ve: view direction
// Input alpha_x, alpha_y: roughness parameters
// Input U1, U2: uniform random numbers
// Output Ne: normal sampled with PDF D_Ve(Ne) = G1(Ve) * max(0, dot(Ve, Ne)) * D(Ne) / Ve.z
vec3 SampleGGXVNDF(vec3 Ve, float alpha_x, float alpha_y, float U1, float U2)
{
    vec3 Vh = normalize(vec3(alpha_x * Ve.x, alpha_y * Ve.y, Ve.z));

    float lensq = Vh.x * Vh.x + Vh.y * Vh.y;
    vec3 T1 = lensq > 0 ? vec3(-Vh.y, Vh.x, 0) * inversesqrt(lensq) : vec3(1,0,0);
    vec3 T2 = cross(Vh, T1);

    float r = sqrt(U1);
    float phi = 2.0 * PI * U2;
    float t1 = r * cos(phi);
    float t2 = r * sin(phi);
    float s = 0.5 * (1.0 + Vh.z);
    t2 = (1.0- s)*sqrt(1.0- t1*t1) + s*t2;

    vec3 Nh = t1*T1 + t2*T2 + sqrt(max(0.0, 1.0- t1*t1- t2*t2))*Vh;

    vec3 Ne = normalize(vec3(alpha_x * Nh.x, alpha_y * Nh.y, max(0.0, Nh.z)));
    return Ne;
}

不过需要注意的是，论文中函数的输入参数 Ve 表示的是局部空间中的视线方向。因此，在调用该函数之前，还需要先完成一次坐标系变换，把世界空间下的视线方向转换到局部空间中：

vec3 WorldToLocal(vec3 v, vec3 T, vec3 B, vec3 N)
{
    return vec3(dot(v, T), dot(v, B), dot(v, N));
}

vec3 LocalToWorld(vec3 v, vec3 T, vec3 B, vec3 N)
{
    return v.x * T + v.y * B + v.z * N;
}

vec3 SampleGGXVNDF_World(vec3 Ve_world, vec3 N, float alpha_x, float alpha_y, float U1, float U2)
{
    N = normalize(N);
    vec3 helper = vec3(1, 0, 0);
    if(abs(N.x)>0.999) helper = vec3(0, 0, 1);
    vec3 T = normalize(cross(N, helper));
    vec3 B = normalize(cross(N, T));

    vec3 Ve_local = WorldToLocal(Ve_world, T, B, N);

    vec3 Vh = normalize(vec3(alpha_x * Ve_local.x, alpha_y * Ve_local.y, Ve_local.z));
    if (Vh.z < 0.0) Ve_local = -Ve_local;

    alpha_x = max(1e-3, alpha_x);
    alpha_y = max(1e-3, alpha_y);
    vec3 Ne_local = SampleGGXVNDF(Ve_local, alpha_x, alpha_y, U1, U2);
    vec3 Ne_world = LocalToWorld(Ne_local, T, B, N);
    return Ne_world;
}

在这些准备工作完成之后，只需要对原有的采样函数和 PDF 评估函数做少量修改，VNDF 采样就可以顺利接入现有框架。为了便于对比和调试，直接加入一个全局开关变量，用来控制是否启用 VNDF 采样：

uniform int uEnableVNDF;

float PdfGlassReflectionIS(vec3 wi, vec3 wo, vec3 Ns, float eta_wo2wi, float alpha)
{
	...

    float pdfH = 1.0;
    if(uEnableVNDF == 1)
        pdfH = smithG1_GGX(NdotV, alpha) * GTR2(NdotH, alpha) * VdotH / NdotV;
    else 
        pdfH = GTR2(NdotH, alpha) * NdotH;

    float Fg = FrDielectric(dot(wo, H), eta_wo2wi);

    return  Fg * pdfH / max(4.0 * VdotH, 1e-6);
}

float PdfGlassTransmissionIS(vec3 wi, vec3 wo, vec3 Ns, float eta_wo2wi, float alpha)
{
	...

    float pdfH = 1.0;
    if(uEnableVNDF == 1)
        pdfH = smithG1_GGX(NdotV, alpha) * GTR2(NdotH, alpha) * VdotH / NdotV;
    else 
        pdfH = GTR2(NdotH, alpha) * NdotH;
    float Fg = FrDielectric(dot(wo, H), eta_wo2wi);

    float sqrtDenom = LdotH + VdotH / eta_wo2wi;
    float dwh_dwi = abs((LdotH) / max(sqrtDenom * sqrtDenom, 1e-6));

    return  (1 - Fg) * pdfH * dwh_dwi;
}

SampleBSDFResult SampleBSDF_Glass(vec3 wo, vec3 Ns, bool bIsInside, int MaterialSlotID)
{
	...
	
    vec3 H;
    if(uEnableVNDF == 1)
        H = SampleGGXVNDF_World(wo, Ns, alpha, alpha, rand(), rand());
    else
        H = SampleGTR2(rand(), rand(), Ns, alpha);

    if (dot(wo, H) < 0.0)
        H = -H;

	...
}

最后补充一点：虽然这里主要展示的是折射分支中的修改方式，但 Disney BRDF 分支中的高光项同样也应该做对应调整。至于清漆项，由于它使用的是 GTR1 法线分布，而不是论文中讨论的 GGX 分布，因此这里并不能直接套用同一套 VNDF 采样方法，暂时保持不变即可。

在附上几张渲染图：

虽然本文已经成功把玻璃材质的 BSDF 整合进了 Disney BRDF，但距离一个更完整的 Disney BSDF 其实还有一步。真正的 Disney BSDF 中，还包含了 Pixar 在 2015 年提出的 BSSRDF 次表面散射模型。

前段时间我也尝试着自己实现了一个版本，不过实际效果只能说渲了依托，而且跑的巨慢，一个正方体就快把我的GPU榨干了。

虽然这里的结果看起来和Blender中次表面着色器Christensen-Burley分支挺接近的，但是视觉效果完全比不上Random Walk分支，思来想去还是用原来的次表面近似吧。

下图是Blender的Christensen-Burley分支

完整代码仓库：
OnikataHoshio/KhRenderer

结语

真是花了好长时间写这些东西啊。

从最开始在 CPU 上写光线追踪，到一步步搭起 GPU 渲染框架；从 CPU 构建 BVH 到 GPU 构建 BVH；再到微表面模型、BRDF、BSDF、体积渲染、BSSRDF、重要性采样、Sobol 序列……一路下来，算是狠狠补了一轮图形学和渲染相关的理论知识。很多一开始看上去有些望而生畏的东西，也是在不断实现和踩坑的过程中，慢慢变得清晰和熟悉起来。

说实话，我也不知道在AI盛行的现在，学习这些东西是否还有意义。也许真的如大多数人所说，这些总有一天会被AI替代吧。不过探索的过程始终是有趣的，如果学习这些内容能感到有趣和满足的话也不失为一种不错的答案。

虽然这些内容大多是几年前甚至更早就已经提出和完善的技术了，但对于现在的我来说，把这一整条线走下来，本身就是一个相当完整的阶段。而作为这一阶段的落幕，这篇文章也算得上一个相当不错的句号了。

其实我之前也想过，要不要顺手把体积渲染相关的内容也写一写，不过现在打消了这个内容的。一方面，用 Ray Marching 去实现基础体积效果，本身并不算特别复杂；另一方面，我也越来越觉得，继续整理这些“复现”和“实现”类的内容，已经很难再带来最开始那种满足感了。或许接下来，比起继续追着已有的内容往前补，尝试去做一些真正属于自己的东西更加能够自娱自乐吧。

参考文献

1. SIGGRAPH 2015 Course: Physically Based Shading in Theory and Practice
2. UCSD CSE 272 Assignment 1:Disney Principled BSDF
3. Sampling the GGX Distribution of Visible Normals
4. Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals
5. Physically Based Rendering: From Theory to Implementation
6. AKGWSB/EzRT: Easy Ray Tracing, a lite renderer and tutorial from theory to implement, with OpenGL
7. DisneyBRDF