Created by potrace 1.15, written by Peter Selinger 2001-2017 Bo's Blog

PBRT读书笔记11-Surface Reflection

 
  • May 15, 2022

这个系列是Physically Based Rendering: From Theory To Implementation的读书笔记，本节主要介绍光线追踪的相关知识。

在前面介绍过的所有辐照度、Monte Carlo积分等知识的基础上，我们终于可以开始实现具体的渲染算法了。回忆一下散射方程的相关概念，基于Monte Carlo积分我们可以对物体表面反(散)射出的radiance进行估计:

\[\begin{aligned} L_o(p, \omega_o) &= \int_{S^2} f (p, \omega_o, \omega_i) L_i(p, \omega_i) \vert \cos{\theta_i} \vert \ d \omega_i \\ &\approx \frac{1}{N} \sum_{j=1}^N \frac{f (p, \omega_o, \omega_j) L_i(p, \omega_j) \vert \cos{\theta_j} \vert}{p(\omega_j)} \end{aligned}\]

在实际采样时我们一般会先从BRDF中采样光线，然后再从场景光源来采样光线。最后将这些光线样本按照一定的权重混合起来，这样的采样方法称为多重重要性采样(multiple importance sampling, MIS)。

Sampling Reflection Functions

首先我们来考虑从BRDF进行采样的情况，在BxDF类中通过Sample_f()成员函数来实现采样的功能。具体来说BxDF::Sample_f()接收出射光线方向wo作为参数，首先会采样出一个入射方向wi，然后利用这两个方向再分别计算光线的概率以及对应的BxDF函数值。在默认实现中，我们使用了cosine weighted hemisphere sampling的方法来从半球上进行采样，光线的概率会存储到pdf中。最后，函数返回入射和出射方向上的BxDF函数值。整个BxDF::Sample_f()函数的定义如下：

Spectrum BxDF::Sample_f(const Vector3f &wo, Vector3f *wi,
        const Point2f &u, Float *pdf, BxDFType *sampledType) const {
    // cosine-sample the hemisphere, flipping the direction if necessary
    *wi = CosineSampleHemisphere(u);
    if (wo.z < 0) wi->z *= -1;

    *pdf = Pdf(wo, *wi);
    return f(wo, *wi);
}

Float BxDF::Pdf(const Vector3f &wo, const Vector3f &wi) const {
    return SameHemisphere(wo, wi) ? AbsCosTheta(wi) * InvPi : 0;
}

inline bool SameHemisphere(const Vector3f &w, const Vector3f &wp) {
    return w.z * wp.z > 0;
}

上面介绍的这种采样方法对于Lambertian BRDF以及Oren–Nayar材质都有很好的表现，对于这些材质模型我们无需进行重载。

Microfacet BxDFs

对于微表面材质，我们可以利用法向分布函数$D(\omega_h)$来进行采样。实际上Torrance–Sparrow BSDF的函数形状基本由几何项$D(\omega_h)$来控制，因此直接从$D(\omega_h)$进行采样是更高效的方法。MicrofacetDistribution基类包含一个Sample_wh()函数用来实现对法向进行采样：

virtual Vector3f Sample_wh(const Vector3f &wo,
                           const Point2f &u) const = 0;

对$D(\omega_h)$进行采样时可以直接对半球面上的所有法向进行采样。以Beckmann–Spizzichino分布为例，利用各向同性的性质可以将Beckmann–Spizzichino分布的PDF进行分解$p_h(\omega_h) = p_h(\theta) \cdot p_h(\phi)$。对于$\phi$，我们直接按照均匀分布进行采样：

\[\phi = 2 \pi \xi\]

其中$\xi$是(0, 1)区间上的一个随机数。接着，根据$p_h(\theta)$进行采样：

\[p_h(\theta) = \frac{2 e^{-\tan^2 \theta / \alpha^2} \sin \theta}{\alpha^2 \cos^3 \theta}\]

这里使用逆变换采样的方法进行处理，首先计算它的CDF：

\[P_j(\theta') = \int_0^{\theta'} \frac{2 e^{-\tan^2 \theta / \alpha^2} \sin \theta}{\alpha^2 \cos^3 \theta} \ d \theta = 1 -e^{-\tan^2 \theta' / \alpha^2}\]

因此我们只需要令随机数$\xi$满足：

\[\xi = 1 -e^{-\tan^2 \theta' / \alpha^2}\]

进而得到：

\[\tan^2 \theta' = -\alpha^2 \ln (1 - \xi)\]

采样出$\phi$和$\tan^2 \theta$后就可以利用三角函数的变换公式来计算出法向$\omega_h$了。

另一种采样方法是考虑微表面的可见性，只对可见的法向进行采样。此时根据masking的理论我们可以重新定义$D(\omega_h)$的分布：

\[D_\omega (\omega_h) = \frac{D(\omega_h) \ G_1(\omega, \omega_h) \ \max (0, \omega \cdot \omega_h)}{\cos \theta}\]

其中$G_1$项表示微表面之间的自遮挡，$\frac{\max (0, \omega \cdot \omega_h)}{\cos \theta}$项则描述了不可见的背面以及投影后的面积。

从修正后的$D_\omega (\omega_h)$进行采样效率更高而且具有更小的方差。具体的采样方法可以参考源码如下：

static Vector3f BeckmannSample(const Vector3f &wi, Float alpha_x, Float alpha_y,
                               Float U1, Float U2) {
    // 1. stretch wi
    Vector3f wiStretched =
        Normalize(Vector3f(alpha_x * wi.x, alpha_y * wi.y, wi.z));

    // 2. simulate P22_{wi}(x_slope, y_slope, 1, 1)
    Float slope_x, slope_y;
    BeckmannSample11(CosTheta(wiStretched), U1, U2, &slope_x, &slope_y);

    // 3. rotate
    Float tmp = CosPhi(wiStretched) * slope_x - SinPhi(wiStretched) * slope_y;
    slope_y = SinPhi(wiStretched) * slope_x + CosPhi(wiStretched) * slope_y;
    slope_x = tmp;

    // 4. unstretch
    slope_x = alpha_x * slope_x;
    slope_y = alpha_y * slope_y;

    // 5. compute normal
    return Normalize(Vector3f(-slope_x, -slope_y, 1.f));
}

static void BeckmannSample11(Float cosThetaI, Float U1, Float U2,
                             Float *slope_x, Float *slope_y) {
    /* Special case (normal incidence) */
    if (cosThetaI > .9999) {
        Float r = std::sqrt(-std::log(1.0f - U1));
        Float sinPhi = std::sin(2 * Pi * U2);
        Float cosPhi = std::cos(2 * Pi * U2);
        *slope_x = r * cosPhi;
        *slope_y = r * sinPhi;
        return;
    }

    /* The original inversion routine from the paper contained
       discontinuities, which causes issues for QMC integration
       and techniques like Kelemen-style MLT. The following code
       performs a numerical inversion with better behavior */
    Float sinThetaI =
        std::sqrt(std::max((Float)0, (Float)1 - cosThetaI * cosThetaI));
    Float tanThetaI = sinThetaI / cosThetaI;
    Float cotThetaI = 1 / tanThetaI;

    /* Search interval -- everything is parameterized
       in the Erf() domain */
    Float a = -1, c = Erf(cotThetaI);
    Float sample_x = std::max(U1, (Float)1e-6f);

    /* Start with a good initial guess */
    // Float b = (1-sample_x) * a + sample_x * c;

    /* We can do better (inverse of an approximation computed in
     * Mathematica) */
    Float thetaI = std::acos(cosThetaI);
    Float fit = 1 + thetaI * (-0.876f + thetaI * (0.4265f - 0.0594f * thetaI));
    Float b = c - (1 + c) * std::pow(1 - sample_x, fit);

    /* Normalization factor for the CDF */
    static const Float SQRT_PI_INV = 1.f / std::sqrt(Pi);
    Float normalization =
        1 /
        (1 + c + SQRT_PI_INV * tanThetaI * std::exp(-cotThetaI * cotThetaI));

    int it = 0;
    while (++it < 10) {
        /* Bisection criterion -- the oddly-looking
           Boolean expression are intentional to check
           for NaNs at little additional cost */
        if (!(b >= a && b <= c)) b = 0.5f * (a + c);

        /* Evaluate the CDF and its derivative
           (i.e. the density function) */
        Float invErf = ErfInv(b);
        Float value =
            normalization *
                (1 + b + SQRT_PI_INV * tanThetaI * std::exp(-invErf * invErf)) -
            sample_x;
        Float derivative = normalization * (1 - invErf * tanThetaI);

        if (std::abs(value) < 1e-5f) break;

        /* Update bisection intervals */
        if (value > 0)
            c = b;
        else
            a = b;

        b -= value / derivative;
    }

    /* Now convert back into a slope value */
    *slope_x = ErfInv(b);

    /* Simulate Y component */
    *slope_y = ErfInv(2.0f * std::max(U2, (Float)1e-6f) - 1.0f);
}

把上面介绍过的两种采样方法整合到一起我们可以得到对Beckmann分布进行采样的函数BeckmannDistribution::Sample_wh()：

Vector3f BeckmannDistribution::Sample_wh(const Vector3f &wo,
        const Point2f &u) const {
    if (!sampleVisibleArea) {
        // Sample full distribution of normals for Beckmann distribution
        // Compute tan2Theta and phi for Beckmann distribution sample
        Float tan2Theta, phi;
        if (alphax == alphay) {
            Float logSample = std::log(1 - u[0]);
            if (std::isinf(logSample)) logSample = 0;
            tan2Theta = -alphax * alphax * logSample;
            phi = u[1] * 2 * Pi;
        } else {
            // Compute tan2Theta and phi for anisotropic Beckmann distribution
            Float logSample = std::log(u[0]);
            phi = std::atan(alphay / alphax *
                            std::tan(2 * Pi * u[1] + 0.5f * Pi));
            if (u[1] > 0.5f)
                phi += Pi;
            Float sinPhi = std::sin(phi), cosPhi = std::cos(phi);
            Float alphax2 = alphax * alphax, alphay2 = alphay * alphay;
            tan2Theta = -logSample /
                         (cosPhi * cosPhi / alphax2 + sinPhi * sinPhi / alphay2);
        }

        // Map sampled Beckmann angles to normal direction wh
        Float cosTheta = 1 / std::sqrt(1 + tan2Theta);
        Float sinTheta = std::sqrt(std::max((Float)0, 1 - cosTheta * cosTheta));
        Vector3f wh = SphericalDirection(sinTheta, cosTheta, phi);
        if (!SameHemisphere(wo, wh)) wh = -wh;

        return wh;
    } else {
        // Sample visible area of normals for Beckmann distribution
        Vector3f wh;
        bool flip = wo.z < 0;
        wh = BeckmannSample(flip ? -wo : wo, alphax, alphay, u[0], u[1]);
        if (flip) wh = -wh;
        return wh;
    }
}

除了采样法向外我们还需要计算法向的pdf，这里需要根据在采样时是否考虑可见性来选择相应的函数：

Float MicrofacetDistribution::Pdf(const Vector3f &wo,
        const Vector3f &wh) const {
    if (sampleVisibleArea)
        return D(wh) * G1(wo) * AbsDot(wo, wh) / AbsCosTheta(wo);
    else
        return D(wh) * AbsCosTheta(wh);
}

采样出法向$\omega_h$后可以利用反射的几何关系来计算入射方向$\omega_i$。但是这里需要注意的是我们目前进行的采样以及计算出的pdf都是针对法向(半程向量)$\omega_h$的，但实际上我们需要的是入射方向$\omega_i$的pdf，因此我们需要进行一些变换把pdf映射到$\omega_i$对应的空间上。

利用反射的几何关系，我们可以推导出入射方向和法向之间的微分立体角满足：

\[\begin{aligned} \frac{d \omega_h}{d \omega_i} &= \frac{\sin \theta_h \ d\theta_h \ d\phi_h}{\sin \theta_i \ d\theta_i \ d\phi_i} \\ &= \frac{\sin \theta_h \ d\theta_h \ d\phi_h}{\sin 2\theta_i \ 2 d\theta_i \ d\phi_i} \\ &= \frac{\sin \theta_h}{4 \cos \theta_h \ \sin \theta_h} \\ &= \frac{1}{4 \cos \theta_h} \\ &= \frac{1}{4 (\omega_i \cdot \omega_o)} \end{aligned}\]

这样就可以得到pdf之间的变换关系：

\[p(\theta) = \frac{p_h (\theta_h)}{4 (\omega_i \cdot \omega_o)}\]

综合这一节介绍过的所有内容，就可以得到重载后的微表面材质采样函数MicrofacetReflection::Sample_f()：

Spectrum MicrofacetReflection::Sample_f(const Vector3f &wo, Vector3f *wi,
        const Point2f &u, Float *pdf, BxDFType *sampledType) const {
    // Sample microfacet orientation  and reflected direction
    Vector3f wh = distribution->Sample_wh(wo, u);
    *wi = Reflect(wo, wh);
    if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);

    // Compute PDF of wi for microfacet reflection
    *pdf = distribution->Pdf(wo, wh) / (4 * Dot(wo, wh));

    return f(wo, *wi);
}

折射的情况与反射类似，我们先采样出法向$\omega_h$然后利用折射的几何关系计算出入射方向$\omega_i$。如果验算发现不会发生全反射则进一步计算光线样本的pdf：

Spectrum MicrofacetTransmission::Sample_f(const Vector3f &wo,
        Vector3f *wi, const Point2f &u, Float *pdf,
        BxDFType *sampledType) const {
    Vector3f wh = distribution->Sample_wh(wo, u);
    Float eta = CosTheta(wo) > 0 ? (etaA / etaB) : (etaB / etaA);
    if (!Refract(wo, (Normal3f)wh, eta, wi))
        return 0;
    *pdf = Pdf(wo, *wi);
    return f(wo, *wi);
}

Float MicrofacetTransmission::Pdf(const Vector3f &wo,
        const Vector3f &wi) const {
    if (SameHemisphere(wo, wi)) return 0;

    // Compute wh from wo and wi for microfacet transmission
    Float eta = CosTheta(wo) > 0 ? (etaB / etaA) : (etaA / etaB);
    Vector3f wh = Normalize(wo + wi * eta);

    // Compute change of variables dwh_dwi for microfacet transmission
    Float sqrtDenom = Dot(wo, wh) + eta * Dot(wi, wh);
    Float dwh_dwi = std::abs((eta * eta * Dot(wi, wh)) / (sqrtDenom * sqrtDenom));

    return distribution->Pdf(wo, wh) * dwh_dwi;
}

FresnelBlend

FresnelBlend类型的BRDF同时包含了漫反射和镜面反射两部分。在对它进行采样时可以使用分层的策略，即当随机数$\xi$小于0.5时按照漫反射材质进行采样，否则按照微表面材质的BRDF进行采样：

Spectrum FresnelBlend::Sample_f(const Vector3f &wo, Vector3f *wi,
        const Point2f &uOrig, Float *pdf, BxDFType *sampledType) const {
    Point2f u = uOrig;
    if (u[0] < .5) {
        u[0] = 2 * u[0];
        // Cosine-sample the hemisphere, flipping the direction if necessary
        *wi = CosineSampleHemisphere(u);
        if (wo.z < 0) wi->z *= -1;

    } else {
        u[0] = 2 * (u[0] - .5f);
        // Sample microfacet orientation wh and reflected direction wi
        Vector3f wh = distribution->Sample_wh(wo, u);
        *wi = Reflect(wo, wh);
        if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);
    }

    *pdf = Pdf(wo, *wi);
    return f(wo, *wi);
}

类似地，在计算样本的pdf时也需要将两种BRDF的pdf加起来并求平均：

Float FresnelBlend::Pdf(const Vector3f &wo, const Vector3f &wi) const {
    if (!SameHemisphere(wo, wi)) return 0;
    Vector3f wh = Normalize(wo + wi);
    Float pdf_wh = distribution->Pdf(wo, wh);
    return .5f * (AbsCosTheta(wi) * InvPi +
                  pdf_wh / (4 * Dot(wo, wh)));
}

Specular Reflection and Transmission

对于包含$\delta$函数的BRDF在进行采样时需要一些额外的处理。回忆$\delta$函数只在$x=0$处有非0值，因此最直接的处理方式是只在唯一可能的方向上进行采样。但是对于需要同时考虑反射和折射的材质，我们需要决定$\delta$函数是返回反射的方向还是折射的方向。在PBRT中使用FresnelSpecular来描述这种理想镜面材质，它利用Fresnel系数来控制对反射还是折射进行采样，FresnelSpecular::Sample_f()函数的返回值分别为光滑镜面反射和折射情况下的BSDF：

Spectrum FresnelSpecular::Sample_f(const Vector3f &wo,
        Vector3f *wi, const Point2f &u, Float *pdf,
        BxDFType *sampledType) const {
    Float F = FrDielectric(CosTheta(wo), etaA, etaB);
    if (u[0] < F) {
        // Compute specular reflection for FresnelSpecular
        // Compute perfect specular reflection direction
        if (sampledType)
            *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION);
        *pdf = F;
        return F * R / AbsCosTheta(*wi);

    } else {
        // Compute specular transmission for FresnelSpecular
        // Figure out which  is incident and which is transmitted
        bool entering = CosTheta(wo) > 0;
        Float etaI = entering ? etaA : etaB;
        Float etaT = entering ? etaB : etaA;

        // Compute ray direction for specular transmission
        if (!Refract(wo, Faceforward(Normal3f(0, 0, 1), wo), etaI / etaT, wi))
            return 0;

        Spectrum ft = T * (1 - F);
        // Account for non-symmetry with transmission to different medium
        if (mode == TransportMode::Radiance)
            ft *= (etaI * etaI) / (etaT * etaT);

        if (sampledType)
            *sampledType = BxDFType(BSDF_SPECULAR | BSDF_TRANSMISSION);
        *pdf = 1 - F;
        return ft / AbsCosTheta(*wi);
    }
}

Fourier BSDF

PBRT中还介绍了Fourier BSDF的采样方法，这里暂时略过。

Application: Estimating Reflectance

到目前为止我们已经完成了PBRT中大部分BxDF的采样函数，这里介绍一个对BxDF进行采样的应用：我们定义材质的定向半球反射率(hemispherical–directional reflectance)为给定出射方向$\omega_o$后BRDF在所有入射方向$\omega_i$上的积分：

\[\rho_{hd} (\omega_o) = \int_{H^2(\mathbf{n})} f_r(\omega_o, \omega_i) \ \vert \cos \theta_i \vert \ d \omega_i\]

显然我们可以使用Monte Carlo积分来估计$\rho_{hd}$，它的估计公式为：

\[\rho_{hd} (\omega_o) \approx \frac{1}{N} \sum_{j=1}^N \frac{f_r(\omega_o, \omega_j) \ \vert \cos \theta_j \vert}{p(\omega_j)}\]

在PBRT中使用BxDF::rho()函数来实现$\rho_{hd}$的估计：

Spectrum BxDF::rho(const Vector3f &w, int nSamples,
        const Point2f *u) const {
    Spectrum r(0.);
    for (int i = 0; i < nSamples; ++i) {
        // Estimate one term of rho
        Vector3f wi;
        Float pdf = 0;
        Spectrum f = Sample_f(w, &wi, u[i], &pdf);
        if (pdf > 0) r += f * AbsCosTheta(wi) / pdf;
    }
    return r / nSamples;
}

类似地，我们同样可以估计半球反射率(hemispherical–hemispherical reflectance)：

\[\begin{aligned} \rho_{hh} &= \frac{1}{\pi} \int_{H^2(\mathbf{n})} \int_{H^2(\mathbf{n})} f_r(\omega', \omega'') \ \vert \cos \theta' \cos \theta'' \vert \ d \omega' \ d \omega'' \\ &= \frac{1}{\pi N} \sum_{j=1}^N \frac{f_r(\omega_j', \omega_j'') \ \vert \cos \theta_j' \cos \theta_j'' \vert}{p(\omega_j') \ p(\omega_j'')} \end{aligned}\] 
Spectrum BxDF::rho(int nSamples, const Point2f *u1,
        const Point2f *u2) const {
    Spectrum r(0.f);
    for (int i = 0; i < nSamples; ++i) {
        // Estimate one term of rho
        Vector3f wo, wi;
        wo = UniformSampleHemisphere(u1[i]);
        Float pdfo = UniformHemispherePdf(), pdfi = 0;
        Spectrum f = Sample_f(wo, &wi, u2[i], &pdfi);
        if (pdfi > 0)
            r += f * AbsCosTheta(wi) * AbsCosTheta(wo) / (pdfo * pdfi);
    }
    return r / (Pi * nSamples);
}

Sampling BSDFs

在BxDF类的基础上我们可以实现BSDF类的采样函数。具体来说，每个BSDF实例包含$N$个独立的BxDF，因此在进行采样时需要分别调用这些BxDF的采样函数然后对它们求平均作为整体的pdf：

\[p(\omega) = \frac{1}{N} \sum_{i=1}^N p_i(\omega)\]

BSDF类的采样函数实现如下：

Spectrum BSDF::Sample_f(const Vector3f &woWorld, Vector3f *wiWorld,
        const Point2f &u, Float *pdf, BxDFType type,
        BxDFType *sampledType) const {
    // Choose which BxDF to sample
    int matchingComps = NumComponents(type);
    if (matchingComps == 0) {
        *pdf = 0;
        return Spectrum(0);
    }
    int comp = std::min((int)std::floor(u[0] * matchingComps), 
                        matchingComps - 1);
    // Get BxDF pointer for chosen component
    BxDF *bxdf = nullptr;
    int count = comp;
    for (int i = 0; i < nBxDFs; ++i)
        if (bxdfs[i]->MatchesFlags(type) && count-- == 0) {
            bxdf = bxdfs[i];
            break;
        }
    
    // Remap BxDF sample u to [0, 1)^2
    Point2f uRemapped(u[0] * matchingComps - comp, u[1]);

    // Sample chosen BxDF
    Vector3f wi, wo = WorldToLocal(woWorld);
    *pdf = 0;
    if (sampledType) *sampledType = bxdf->type;
    Spectrum f = bxdf->Sample_f(wo, &wi, uRemapped, pdf, sampledType);
    if (*pdf == 0) 
        return 0;
    *wiWorld = LocalToWorld(wi);

    // Compute overall PDF with all matching BxDFs
    if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1)
        for (int i = 0; i < nBxDFs; ++i)
            if (bxdfs[i] != bxdf && bxdfs[i]->MatchesFlags(type))
                *pdf += bxdfs[i]->Pdf(wo, wi);
    if (matchingComps > 1) *pdf /= matchingComps;

    // Compute value of BSDF for sampled direction
    if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1) {
        bool reflect = Dot(*wiWorld, ng) * Dot(woWorld, ng) > 0;
        f = 0.;
        for (int i = 0; i < nBxDFs; ++i)
            if (bxdfs[i]->MatchesFlags(type) &&
                ((reflect && (bxdfs[i]->type & BSDF_REFLECTION)) ||
                 (!reflect && (bxdfs[i]->type & BSDF_TRANSMISSION))))
                f += bxdfs[i]->f(wo, wi);
    }
    return f;
}

BSDF::Pdf()函数与采样函数类似，同样需要对自身包含的BxDF进行遍历然后再求平均。

Sampling Light Sources

Lights with Singularities

Sampling Shapes

Area Lights

Infinite Area Lights

Direct Lighting

Integrator Interface and SamplerIntegrator

Estimating the Direct Lighting Integral

The Light Transport Equation

Basic Derivation

Analytic Solutions to the LTE

The Surface Form of the LTE

Integral over Paths

Delta Distributions in the Integrand

Partitioning the Integrand

Path Tracing

Overview

Path Sampling

Incremental Path Construction

Implementation

Reference

Search