PBRT笔记(13)——光线传播1：表面反射

采样反射函数

BxDF::Sample_f()方法根据与相应的散射函数相似的分布来选择方向。在8.2节中，该方法用于寻找来自完美镜面的反射和透射光线;在这里讲介绍实现其他类型的采样技术。

BxDF::Sample_f()在[0,1)范围内取得两个使用反演法取得的样本，其中这些样本是使用分层采样或者低偏差采样生成的，所以这些样本本身具有很好的分布性。

这个方法默认采用余弦加权的半球采样。这个样本分布对于大多数BRDF会产生正确的结果。

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);
}

CosineSampleHemisphere()函数返回一个在半球中绕着反射坐标系中(0,0,1)向量的方向向量。如果ωo在另一种半球中，那么ωi就会被翻转到与ωo相同的半球中。所以pbrt中无需手动翻转不同半球的方向向量向量（因为这里翻转了）。

注意：如果重新实现了BxDF::Sample_f()，那BxDF::Sample_f()也必须重新实现，以便返回一样的pdf。默认的Pdf()会先去检测入射方向与出射方向是否在同一半球，如果不在，则返回0。

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与OrenNayar模型则无需重新实现这两个函数。

微表面BxDFs

8.4节中定义的以微表面为基础的反射模型是基于一个分布微表面几何分布函数 D(ωh),其中每个微表面表现出完美的镜面反射或者透射。

因为D(ωh)函数主要负责Torrance-Sparrow BSDF中微表面几何分布描述,所以使用基于抽样的方法来描述其分布是相当有效的。

该方法首先从微面分布中抽样，从而取得特定微面的方向，然后利用镜面反射或透射公式求出入射方向。经典的微表面方向采样方法是直接采样D(ωh)。

MicrofacetDistribution实现了以上想法，并且它存储了一个布尔值用于决定是否使用采样技术。因为在实际运用中，基于可见微表面采样方法比基于整体采样的方法更有效。

Beckmann–Spizzichino分布的所有法线分布的采样方法返回球面坐标系的tan2θ与φ，并转化成标准化的方向向量wh。为了采样，我们将推导球面坐标系中下的表达式。因为是各项同性分布，它是独立于φ的，所以这个分布的pdf ph(θ,φ)可以分离成ph(θ)和ph(φ)

ph(φ)为常数1/2π，所以φ可以被采样为：

$\phi=2\pi\xi$

对于p(θh)：

$p(\theta_h)=\frac{e^{-tan^2 \theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}$

知乎有人说这里pbrt错了,应为: $D(\omega_h)=\frac{e^{-tan^2 \theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}$ 链接：https://zhuanlan.zhihu.com/p/64080095

其中α是粗糙度系数。我们可以使用反演方法，对给定的均匀随机ξ，从这个分布中采样方向θ。首先从公式中获取pdf，之后再积分得到cdf：

$P(\theta')=\int^{\theta'}_{0}\frac{e^{-tan^2 \theta_h/\alpha^2}}{\pi\alpha^2cos^4\theta_h}d\theta_h=1-e^{-tan^2 \theta_h'/\alpha^2}$

为了采样，我们需要解以下式子：

$\xi=1-e^{-tan^2 \theta_h'/\alpha^2}$

得：

$tan^2\theta'=-\alpha^2log \xi$

Vector3f BeckmannDistribution::Sample_wh(const Vector3f &wo,const Point2f &u) const {
    if (!sampleVisibleArea) {
        // Sample full distribution of normals for Beckmann distribution

        // Compute $\tan^2 \theta$ and $\phi$ for Beckmann distribution sample
        Float tan2Theta, phi;
        if (alphax == alphay) {
            Float logSample = std::log(1 - u[0]);
            DCHECK(!std::isinf(logSample));
            tan2Theta = -alphax * alphax * logSample;
            phi = u[1] * 2 * Pi;
        } else {
            // Compute _tan2Theta_ and _phi_ for anisotropic Beckmann
            // distribution
            Float logSample = std::log(1 - u[0]);
            DCHECK(!std::isinf(logSample));
            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;
    }
}

实际上基于可见微表面采样方法与基于整体采样的方法的结果相当接近。而基于整体采样的方法的效率不佳。同时结果表明，样本可以直接从方程定义的分布中提取，因为这个分布更好地匹配完整的Torrance-Sparrow模型。(方程(8.19))比D(ωh),有更少的方差结果。

基于从微表面方向分布中采样，实现为MicrofacetReflection::Sample_f()。（之后pbrt推了一下基于微表面的brdf的公式，这里略过）

Spectrum MicrofacetReflection::Sample_f(const Vector3f &wo, Vector3f *wi,const Point2f &u, Float *pdf,BxDFType *sampledType) const {
    // Sample microfacet orientation $\wh$ and reflected direction $\wi$
    if (wo.z == 0) return 0.;
    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);
}

关于透射：

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);
}

菲尼尔融合

FresnelBlend类定义为diffuse与glossy的混合结果。一种相对直接的BRDF采样方法是根据余弦加权与微表面分布进行采样。在计算前先要判断ξ1是否大于或小于0.5，之后再将它重新映射ξ1至[0,1)再做后续处理。

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] = std::min(2 * u[0], OneMinusEpsilon);
        // Cosine-sample the hemisphere, flipping the direction if necessary
        *wi = CosineSampleHemisphere(u);
        if (wo.z < 0) wi->z *= -1;
    } else {
        u[0] = std::min(2 * (u[0] - .5f), OneMinusEpsilon);
        // 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);
}

高光反射与透射

之前使用狄拉克分布函数定义BRDF的高光反射与BTDF的高光透射。尽管如此，当使用这一类采样机制和pdf函数时仍需留意一些规则。狄拉克分布函数的定义为：
$\delta(x)=0$ 对于所有x皆不为0，同时：

$\int^\infty_{-\infty}\delta(x)dx=1$

因此上式可定义为pdf，对于所有x=0含有一个0值。从这样一个分布中生成样本比较简单，因为它只需使用一个可能值。以这种方式思考时，实现SpecularReflection与SpecularTransmission的BxDFs的Sample_f()就能很自然地使用蒙特卡洛采样框架。

确定pdf的返回值并不是想象中的那么简单，严格地说，狄拉克分布函数并不是一个真正的函数，且需要定义为另一个函数的极限——例如：一个宽度趋近于0的单位面积的盒子。据此,δ(0)的值趋近于无穷大。pdf返回一个无穷大或非常大的值会导致渲染器产生错误的结果。

回忆一下，使用含有狄拉克项的BSDFs中的tr()一样含有狄拉克成分。但我们在它们的Sample_f()中忽略了这个细节。因此，此类带有BSDF项的散射方程的蒙特卡罗估计为：

$\frac{1}{N}\sum^N_i\frac{f_r(p,\omega_o,\omega_i)L_i(p,\omega_i)\left|cos\theta_i\right|}{p(\omega_i)}=\frac{1}{N}\sum^N_i\frac{\rho_{hd}(\omega_o)\frac{\phi(\omega-\omega_i)}{\left|cos\theta_i\right|}L_i(p,\omega_i)\left|cos\theta_i\right|}{p(\omega_i)}$

其中ρhd(ωo)返回半球方向反射率，ω是完美镜面反射或者投射的方向。

因为pdf p(ωi)也含有狄拉克项，且 $p(\omega_i)=\delta(\omega-\omega_i)$ ,两个狄拉克项可以抵消，因此估算函数为： $\rho_hd(\omega_o)L_i(p,\omega)$

因此，在使用Sample_f()进行采样镜面反射或者透射时，pdf返回一个常量1。总结一下：当镜面反射BxDFs的pdf中隐含有一个狄拉克项时可以通过BSDF中隐含狄拉克函项进行抵消。因此相应的pdf方法返回0。

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
        *wi = Vector3f(-wo.x, -wo.y, wo.z);
        if (sampledType)
            *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION);
        *pdf = F;
        return F * R / AbsCosTheta(*wi);
    } else {
        // Compute specular transmission for _FresnelSpecular_

        // Figure out which $\eta$ 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);
    }
}

应用：估算反射

首先回忆一下半球反射率：

$\rho_hd(\omega_o)=\int_{H^2(n)}f_r(\omega_o,\omega_i)\left|cos\theta_i\right|d\omega_i$

BxDF::rho()的实现采用了两个额外参数：nSamples（采样数）与u（样本数组，提供样本）。这个方法对任意BxDF计算其蒙特卡洛估算值，使用BxDF的重要性采样进行估算。

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_\roman{hd}$
        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;
}

实际上，估计量的计算就是对反射函数的分布进行采样，之后再除以pdf。

$\frac{1}{N}\sum^N_j\frac{f_r(\omega,\omega_j)\left|cos\theta_j\right|}{p(\omega_j)}$

BxDF::Sample_f()返回ωj、p(ωj)、fr(ωoωj)的值。不过需要注意p(ωj)=0的情况，因为除以0，r会变得无穷大。

半球-半球反射率可以使用类似的方法估算:

$\rho_hh=\frac{1}{\pi}\int_{H^2(n)}\int_{H^2(n)}f_r(\omega',\omega'')\left|cos\theta'cos\theta''\right|d\omega'd\omega''$ $\frac{1}{\pi N}\sum^N_j\frac{f_r(\omega_j',\omega_j'')\left|cos\theta_j'cos\theta_j''\right|}{p(\omega_j')p(\omega_j'')}$

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_\roman{hh}$
        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);
}

BSDFs采样

对这些单个BxDFs采样，我们定义了BSDF类，并专门实现了Sample_f()，该方法由积分器调用，并根据BSDF分布进行采样。BSDF类存储一个或多个BxDFs指针，虽然这些指针对象可以单独采样，但是这里我们将使用它们密度函数的均值进行采样。

$p(\omega)=\frac{1}{N}\sum^N_j P_i(\omega)$

该方法首先确定当前样本该使用哪种BxDF进行采样，这个根据BxDFType进行匹配。

Spectrum BSDF::Sample_f(const Vector3f &woWorld, Vector3f *wiWorld,const Point2f &u, Float *pdf, BxDFType type,BxDFType *sampledType) const {
    ProfilePhase pp(Prof::BSDFSampling);
    //取得匹配的BxDF的数量，NumComponents进行了匹配判断
    int matchingComps =NumComponents(type);
    if (matchingComps == 0) {
        *pdf = 0;
        if (sampledType) *sampledType = BxDFType(0);
        return Spectrum(0);
    }
    int comp =
        std::min((int)std::floor(u[0] * matchingComps), matchingComps - 1);

    //取得匹配的bxdf指针
    BxDF *bxdf = nullptr;
    int count = comp;
    for (int i = 0; i < nBxDFs; ++i)
        if (bxdfs[i]->MatchesFlags(type) && count-- == 0) {
            bxdf = bxdfs[i];
            break;
        }
    CHECK(bxdf != nullptr);
    VLOG(2) << "BSDF::Sample_f chose comp = " << comp << " / matching = " <<
        matchingComps << ", bxdf: " << bxdf->ToString();

    //将样本重新映射至[0,1)区间。
    因为要通过u[0]的样本来确定要使用哪个BxDF组件进行采样（对应的采样函数里判断，比如圆盘采样），所以我们不能要用组件的Sample_f()方法，因为他不再是均匀分布的。
    Point2f uRemapped(std::min(u[0] * matchingComps - comp, OneMinusEpsilon),u[1]);

    //采样对应BxDF，对应BxDF返回对应BxDF坐标系的向量,但最终都需要转换成世界坐标
    Vector3f wi, wo = WorldToLocal(woWorld);
    if (wo.z == 0) return 0.;
    *pdf = 0;
    if (sampledType) *sampledType = bxdf->type;
    Spectrum f = bxdf->Sample_f(wo, &wi, uRemapped, pdf, sampledType);
    VLOG(2) << "For wo = " << wo << ", sampled f = " << f << ", pdf = "
            << *pdf << ", ratio = " << ((*pdf > 0) ? (f / *pdf) : Spectrum(0.))
            << ", wi = " << wi;
    if (*pdf == 0) {
        if (sampledType) *sampledType = BxDFType(0);
        return 0;
    }
    *wiWorld = LocalToWorld(wi);

     //为了计算被采样的ωi的pdf，我们需要对所有所有的BxDF的pdf求均值。（跳过镜面反射类型，因为pdf公式中有狄拉克项，将其他类型的pdf加入是不正确的）
    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;

    //对于给定的采样方向，这个方法需要计算对应的BSDF值（调用对应BxDF::f()），对于镜面反射组件则跳过计算（其BxDF::f()返回0）
    if (!(bxdf->type & BSDF_SPECULAR)) {
        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);
    }
    VLOG(2) << "Overall f = " << f << ", pdf = " << *pdf << ", ratio = "
            << ((*pdf > 0) ? (f / *pdf) : Spectrum(0.));
    return f;
}

对光源进行采样

因为根据辐射度理论，对某一位置根据BSDF分布计算光线照度其效率非常低，所以更好的方法就是直接采样基于光源的分布，其基础方法为

virtual Spectrum Sample_Li(const Interaction &ref, const Point2f &u,
Vector3f *wi, Float *pdf, VisibilityTester *vis) const = 0;

参数u代表了对光源采样后得到的2d样本，选择采样方向的pdf返回为* PDF。

点光源

对于当前镜面反射与透射的采样框架来说，根据狄拉克分布定义的光源可以自然融入。但仍需谨慎处理其中若干采样方法。

该光源通过狄拉克分布来描述，接受点将仅从一个方向接受照明。因此采样问题将变得简单。也因此PointLight::Pdf_Li()将返回0。

在之前的章节中已经贴过代码了（包括了包括了聚光灯、投射光源、方向光等），所以这里就不贴了。

形状采样

在pbrt中，面光源是通过在Shape类上附加发射轮廓来定义的。因此，为了对这些光源发出入射光线进行采样，在Shape的表面上生成样本是很有用的。

有两种形状采样方法，都名为shape::Sample()。

第一种选择Shape表面遵照表面积分布生成采样点进行运算。最终返回有关采样点含有几何信息的Interaction对象。了初始化采样点的位置p和法线n外，Sample()方法应该将Interaction::pError设置为计算的p值中浮点四舍五入误差的边界。其pdf为面积的倒数。

virtual Interaction Sample(const Point2f &u) const = 0;

第二种是将一个被积分的Shape表面的点作为一个参数，这种方法对于照明特别有用：因为调用者可以传入要被点亮的点，并允许Shape实现可见性判断逻辑：
确保他们只对该点可见的部分Shape采样。默认情况下调用第一种方法。

virtual Interaction Sample(const Interaction &ref,const Point2f &u) const {
    return Sample(u);
}

与第一种Shape采样方法不同：

该方法是根据Shape表面积的概率密度概率来生成采样点的。
使用参考点ref的立体来计算角密度概率函数。（对应的pdf方法也与第一种不同）
这种差异来源于面光源采样估算直接照明积分是一个参考点方向的积分。

Float Shape::Pdf(const Interaction &ref, const Vector3f &wi) const {
    // Intersect sample ray with area light geometry
    Ray ray = ref.SpawnRay(wi);
    Float tHit;
    SurfaceInteraction isectLight;
    // Ignore any alpha textures used for trimming the shape when performing
    // this intersection. Hack for the "San Miguel" scene, where this is used
    // to make an invisible area light.
    if (!Intersect(ray, &tHit, &isectLight, false)) return 0;

    // Convert light sample weight to solid angle measure
    Float pdf = DistanceSquared(ref.p, isectLight.p) /
                (AbsDot(isectLight.n, -wi) * Area());
    if (std::isinf(pdf)) pdf = 0.f;
    return pdf;
}

给定一个参考点与方向ωi,Pdf()会以此生成一条射线用于求交测试（只对面光源Shape进行求交测试），如果射线不与Shape相交，pdf会被认为是0。

要计算参考点立面角的pdf的值，首先需要计算基于Shape表面积的pdf，之后再除以因子将其转换：

$\frac{d\omega_i}{dA}=\frac{cos\theta_o}{r^2}$

Float pdf = DistanceSquared(ref.p, isectLight.p) /(AbsDot(isectLight.n, -wi) * Area());

圆盘采样

圆盘采样是使用同心圆盘采样函数是在单位圆盘中找一个采样点，再进行缩放与偏移，使其复合当前圆盘。注意：这个方法不适合局部圆盘（1、内径不为0 2、角度小于2π）。

Interaction Disk::Sample(const Point2f &u) const {
    Point2f pd = ConcentricSampleDisk(u);
    Point3f pObj(pd.x * radius, pd.y * radius, height);
    Interaction it;
    it.n = Normalize((*ObjectToWorld)(Normal3f(0, 0, 1)));
    if (reverseOrientation) it.n *= -1;
    it.p = (*ObjectToWorld)(pObj, Vector3f(0, 0, 0), &it.pError);
    return it;
}

三角形采样

之前章节有介绍UniformSampleTriangle()，返回三角形中基于重心坐标系的均匀的采样点。

Interaction Triangle::Sample(const Point2f &u, Float *pdf) const {
    Point2f b = UniformSampleTriangle(u);
    // Get triangle vertices in _p0_, _p1_, and _p2_
    const Point3f &p0 = mesh->p[v[0]];
    const Point3f &p1 = mesh->p[v[1]];
    const Point3f &p2 = mesh->p[v[2]];
    Interaction it;
    it.p = b[0] * p0 + b[1] * p1 + (1 - b[0] - b[1]) * p2;
    // Compute surface normal for sampled point on triangle
    it.n = Normalize(Normal3f(Cross(p1 - p0, p2 - p0)));
    // Ensure correct orientation of the geometric normal; follow the same
    // approach as was used in Triangle::Intersect().
    if (mesh->n) {
        Normal3f ns(b[0] * mesh->n[v[0]] + b[1] * mesh->n[v[1]] +
                    (1 - b[0] - b[1]) * mesh->n[v[2]]);
        it.n = Faceforward(it.n, ns);
    } else if (reverseOrientation ^ transformSwapsHandedness)
        it.n *= -1;

    // Compute error bounds for sampled point on triangle
    Point3f pAbsSum =
        Abs(b[0] * p0) + Abs(b[1] * p1) + Abs((1 - b[0] - b[1]) * p2);
    it.pError = gamma(6) * Vector3f(pAbsSum.x, pAbsSum.y, pAbsSum.z);
    *pdf = 1 / Area();
    return it;
}

球体采样

与圆盘采样一样，这里介绍的方法不适用于非完整球体。sphere::Sample()仅仅使用UniformSampleSphere()函数在单位球面上生成一个点，并将该点按球面半径进行缩放。

Interaction Sphere::Sample(const Point2f &u, Float *pdf) const {
    Point3f pObj = Point3f(0, 0, 0) + radius * UniformSampleSphere(u);
    Interaction it;
    it.n = Normalize((*ObjectToWorld)(Normal3f(pObj.x, pObj.y, pObj.z)));
    if (reverseOrientation) it.n *= -1;
    // Reproject _pObj_ to sphere surface and compute _pObjError_
    pObj *= radius / Distance(pObj, Point3f(0, 0, 0));
    Vector3f pObjError = gamma(5) * Abs((Vector3f)pObj);
    it.p = (*ObjectToWorld)(pObj, pObjError, &it.pError);
    *pdf = 1 / Area();
    return it;
}

对于给定一个需要被照亮的点，还有另一种球体采样函数。虽然之前的那种对整个球体表面的均匀采样是完全正确的，但是更好的方法是对可见区域的球面点进行采样。这里的采样程序从参考点开始，均匀地采样球面所对的立体角上的方向，然后计算球面上与采样方向相对应的点。

Interaction Sphere::Sample(const Interaction &ref, const Point2f &u,Float *pdf) const {
    Point3f pCenter = (*ObjectToWorld)(Point3f(0, 0, 0));

    //如果在球体内部就进行均匀采样
    //因为从内部观察，整个球体都是可见的，这里需要注意使用的参考点pOrigin是使用OffsetRayOrigin()计算的，这样可以保证参考点来自与球面相交的射线，从而保证点位于球面正确的一面
    Point3f pOrigin =
        OffsetRayOrigin(ref.p, ref.pError, ref.n, pCenter - ref.p);
    if (DistanceSquared(pOrigin, pCenter) <= radius * radius) {
        Interaction intr = Sample(u, pdf);
        Vector3f wi = intr.p - ref.p;
        if (wi.LengthSquared() == 0)
            *pdf = 0;
        else {
            // Convert from area measure returned by Sample() call above to
            // solid angle measure.
            wi = Normalize(wi);
            *pdf *= DistanceSquared(ref.p, intr.p) / AbsDot(intr.n, -wi);
        }
        if (std::isinf(*pdf)) *pdf = 0.f;
        return intr;
    }

    //计算坐标系用于球体采样
    Vector3f wc = Normalize(pCenter - ref.p);
    Vector3f wcX, wcY;
    CoordinateSystem(wc, &wcX, &wcY);

    // Sample sphere uniformly inside subtended cone

    // Compute $\theta$ and $\phi$ values for sample in cone
    Float sinThetaMax2 = radius * radius / DistanceSquared(ref.p, pCenter);
    Float cosThetaMax = std::sqrt(std::max((Float)0, 1 - sinThetaMax2));
    Float cosTheta = (1 - u[0]) + u[0] * cosThetaMax;
    Float sinTheta = std::sqrt(std::max((Float)0, 1 - cosTheta * cosTheta));
    Float phi = u[1] * 2 * Pi;

    // Compute angle $\alpha$ from center of sphere to sampled point on surface
    Float dc = Distance(ref.p, pCenter);
    Float ds = dc * cosTheta -
               std::sqrt(std::max(
                   (Float)0, radius * radius - dc * dc * sinTheta * sinTheta));
    Float cosAlpha = (dc * dc + radius * radius - ds * ds) / (2 * dc * radius);
    Float sinAlpha = std::sqrt(std::max((Float)0, 1 - cosAlpha * cosAlpha));

    // Compute surface normal and sampled point on sphere
    Vector3f nWorld =
        SphericalDirection(sinAlpha, cosAlpha, phi, -wcX, -wcY, -wc);
    Point3f pWorld = pCenter + radius * Point3f(nWorld.x, nWorld.y, nWorld.z);

    // Return _Interaction_ for sampled point on sphere
    Interaction it;
    it.p = pWorld;
    it.pError = gamma(5) * Abs((Vector3f)pWorld);
    it.n = Normal3f(nWorld);
    if (reverseOrientation) it.n *= -1;

    // Uniform cone PDF.
    *pdf = 1 / (2 * Pi * (1 - cosThetaMax));

    return it;
}

如果参考点在球体外面，那么被遮蔽的p与球体所成的角度为：

$\theta_{max}=arcsin(\frac{r}{\left|p-p_c\right|})=arccos \sqrt {1-(\frac{r}{\left|p-p_c\right|})}$

其中r为球体半径，pc为圆心。采样方法使用上面的公式，计算面向球体的余弦角θmax，然后在方向锥内对方向进行均匀采样。

给定一个基于之前计算的采样坐标系的样本角(θ,φ)，我们可以直接计算在球面上对应的样本。这个方法有3个步骤：

如果我们用来表示参考点到球体中心的距离的变量dc与参考点组成角度为θ的直角三角形，那么我们可以计算得出三角形另外两边的长度。
接下来，考虑图14.9(b)所示的直角三角形，其中斜边是长度等于球面半径r。那这个三角形的第三条边为：那么从参考点指向球体采样点的长度为：我们现在可以使用余弦定理计算球体中心与球面采样点采样点所组成的边与图中dc的夹角α
面光源
有了以上的介绍的Shape采样方法，实现DiffuseAreaLight::Sample_Li()就非常简单了。大部分的工作内容都集中在Shape部分，DiffuseAreaLight只需要复制适当的值来输出参数并计算发射的辐射值。
```
Spectrum DiffuseAreaLight::Sample_Li(const Interaction &ref, const Point2f &u,Vector3f *wi, Float *pdf,VisibilityTester *vis) const {
 ProfilePhase _(Prof::LightSample);
 Interaction pShape = shape->Sample(ref, u, pdf);
 pShape.mediumInterface = mediumInterface;
 if (*pdf == 0 || (pShape.p - ref.p).LengthSquared() == 0) {
     *pdf = 0;
     return 0.f;
 }
 *wi = Normalize(pShape.p - ref.p);
 *vis = VisibilityTester(ref, pShape);
 return L(pShape, -*wi);
}
```
相对的Pdf_Li()直接返回Shape::Pdf()的值得以实现。
无限大的面光源
使用带有环境全景贴图的InfiniteAreaLights通常在不同的方向上有着较大的亮度差异：例如，白天的天空环境全景图。其差异会大到几千倍的差距，通过对InfiniteAreaLights进行采样并匹配照度分布，可以充分地减少图像的方差，减少最终结果的噪点。主要的采样步骤为：
我们在图片坐标系(u,v)下，定义一个二维的分段-常量分布函数p(u,v)，该分布函数对应了环境全景图的辐射度部分。
我们将采用13.6.7节中介绍的采样方法将二位样本转化成从p(u,v)分布中抽取的样本。
我们根据概率密度除以(u, v)定义了单位球方向上的概率密度函数。

第一步是使用Spectrum::y()将使用环境全景图像素定义的连续不断的光谱辐射度函数转换为分段常量函数，之后再来计算采样点的亮度值。在计算过程中有3点需要注意：

这个函数计算与原始环境贴图中像素数相同数量相同的点辐射函数值，但是，它可以使用更多或更少的采样点，这也会相应地增多或减少内存的使用，但仍然生成有效的抽样分布。但是如果抽样分布与函数不匹配，那么将会浪费内存，且增量收益很小。
分段常量函数值被存储在图片中，是通过找到MIPMap::Lookup()（不仅仅是复制对应的值）轻微模糊辐射度函数后找到的值。图14.11展示了这种想法。由于我们使用的是分段常数函数而不是分段线性函数进行采样，因此它必须考虑到这个问题，以确保在任何辐射度函数非零的点进行采样的概率大于零。

每个在缓存中的图片像素值都会被乘以对应的的sinθ值（对应球体经纬度映射的θ）。注意：这种乘法对于抽样方法的正确性没有影响。因为sinθ的值总是大于0(0,π)范围,我们只是重塑抽样的PDF。这里调整PDF是为了消除在采样方法中使用二维图像映射到单位球面所造成的失真影响;

 //这段在构造函数中
 // Compute scalar-valued image _img_ from environment map
 int width = 2 * Lmap->Width(), height = 2 * Lmap->Height();
 std::unique_ptr<Float[]> img(new Float[width * height]);
 float fwidth = 0.5f / std::min(width, height);
 ParallelFor(
     [&](int64_t v) {
         Float vp = (v + .5f) / (Float)height;
         Float sinTheta = std::sin(Pi * (v + .5f) / height);
         for (int u = 0; u < width; ++u) {
             Float up = (u + .5f) / (Float)width;
             img[u + v * width] = Lmap->Lookup(Point2f(up, vp), fwidth).y();
             img[u + v * width] *= sinTheta;
         }
     },
     height, 32);
 // Compute sampling distributions for rows and columns of image
 distribution.reset(new Distribution2D(img.get(), width, height));

有这些预计算数据，抽样部分的任务就比较简单了。给定一个样本u /[0,1)2，它使用第13.6.7节中描述的抽样算法从函数的分布中抽取一个样本，该算法给出一个(u,v)值和获取该样本的PDF的值p(u, v)。

Spectrum InfiniteAreaLight::Sample_Li(const Interaction &ref, const Point2f &u,Vector3f *wi, Float *pdf,VisibilityTester *vis) const {
 ProfilePhase _(Prof::LightSample);
 // Find $(u,v)$ sample coordinates in infinite light texture
 Float mapPdf;
 Point2f uv = distribution->SampleContinuous(u, &mapPdf);
 if (mapPdf == 0) return Spectrum(0.f);

 // Convert infinite light sample point to direction
 Float theta = uv[1] * Pi, phi = uv[0] * 2 * Pi;
 Float cosTheta = std::cos(theta), sinTheta = std::sin(theta);
 Float sinPhi = std::sin(phi), cosPhi = std::cos(phi);
 *wi =
     LightToWorld(Vector3f(sinTheta * cosPhi, sinTheta * sinPhi, cosTheta));

 // Compute PDF for sampled infinite light direction
 *pdf = mapPdf / (2 * Pi * Pi * sinTheta);
 if (sinTheta == 0) *pdf = 0;

 // Return radiance value for infinite light direction
 *vis = VisibilityTester(ref, Interaction(ref.p + *wi * (2 * worldRadius),
                                          ref.time, mediumInterface));
 return Spectrum(Lmap->Lookup(uv), SpectrumType::Illuminant);
}

(u,v)样本映射至球体坐标：

$(\theta,\phi)=(\pi v,2\pi u)$

回顾一下，光源采样程序返回概率密度值，必须根据在单位球面上的测量立面角来定义。因此，这个程序必须马上计算所使用处于 $[0,2]^2$ 的图像函数采样密度与使用经纬度映射将图像映射到单位球面后的相应密度之间的转换。

雅可比矩阵的行列式的绝对值|Jg|=2π2。应用13.5.1的多维变化方程，我们可以得到球体密度坐标(θ,φ)。

$p(\theta,\phi)=\frac{p(u,v)}{2\pi^2}$

从球体坐标可以很容易确定从(r,θ,φ)映射到(x, y,z)的雅可比矩阵的行列式的绝对值是 $r^2sin\theta$ 。因为我们对单位球感兴趣，r=1，所以：

$p(\omega)=\frac{p(\theta,\phi)}{sin\theta}=\frac{p(u,v)}{2\pi^2sin\theta}$

应用了这项技术的主要原因是：它让我们从图像映射定义的分段常数分布中采样，并将样本及其概率密度转换为单位球面上的方向。

我们现在可以看到为什么初始化程序将分段常量采样函数乘以sinθ。例如：常量环境全景图使用p(u,v)采样技术，所有的(θ,φ)被选中的可能是相同的。然而，由于是在球体上映射方向，这将导致球面极点附近的采样过多，而不是对球面所有方向进行平均采样来得到一个更理想的结果。1/sinθ的pdf将纠正这种非均匀的采样方向，以便在蒙特卡洛估算中得到正确的结果。鉴于这种情况，最好修改一下p(u, v)抽样分布，这样减少了采样时选择到球体极点的可能。

直接照明

之前我们介绍了光线传输方程的一般性质，现在我们将实现DirectLightingIntegrator类（直接照明积分器），只包含光线从光源传输到着色点的直接光照，忽略非发光物的间接光照（但不包括基础的镜面反射与透射效果）。

该类提供了两种不同策略来计算直接照明。每种方法都会通过计算给出一个该点坐标给定方向向外辐射的无偏差的估算结果。通过LightStrategy枚举来选择方法。

第一种UniformSampleAll：遍历所有的灯光，并从灯光处取得基于Light::nSamples量的采样样本，并对结果进行求和。
第二种UniformSampleOne：只从一个灯光中随机取一个样本。

根据渲染场景的具体情况，这两种方法都可能合适。例如每个像素上拥有许多图片样本（例如：解决景深没有过多的噪点）那么UniformSampleOne将更为合适：每像素上所有图片样本将会采样直接照明从而得到足够高质量的结果。但是如果每像素的样本太少那么最好对所有的灯光进行采样，以获得无噪点的结果。(PS.这里的意思是如果你MSAA高，灯光采样数就可以降低，但是结果还是一样的)

DirectLightingIntegrator的构造函数会传递摄像机与采样积分器，以及最大光线递归数与直接照明策略。积分器所需样本的数量与类型取决于使用的采样策略：如果一个从单个光源中取得的单独样本，然后通过Sampler::Get2D()获取两组足够多的二维样本，其中一个用于记录光源样本（选择光源上的坐标），另一个用于从BSDF采样散射方向。

当从每个灯光取得多个样本时，积分器会在主进程处理前请求样本数组。使用样本数组能更好地分别调用Sampler::Get2D()。因为它是的采样器有机会使用改进的样本放置技术，即在数组中优化样本的分布。从而对在着色点所有光源的样本进行优化。

void DirectLightingIntegrator::Preprocess(const Scene &scene,
Sampler &sampler) {
if (strategy == LightStrategy::UniformSampleAll) {
<Compute number of samples to use for each light 853>
<Request samples for sampling all lights 853>
}
}

对于LightStrategy::UniformSampleAll策略，每个光源的nSamples成员变量存储了需要的采样数。但是这里的积分器只是用这个值作为起点：Sampler::RoundCount()会通过特定的样本生成技术将其更改为更加合适的值。每个灯光最终的样本数将记录在nLightSamples成员变量。

现在样本数组可以被请求了，这里有两个重要细节：

尽管对所有射线最大深度范围内都提出了单独的样本请求，但这并不意味着所有相交点都会获得样本数组。相反，一旦调用Sampler::Get2DArray()检测到样本数组，该数组将不再返回。如果存在镜面反射与透射，那么可能有 $2^{maxDepth+1}-1$ 个与摄像机射线相交的相交点。如果样本数组耗尽，积分器将会切换成UniformSampleOne模式。
注意请求数组顺序，他们将被积分器使用：每个相交点，2个数组会被用于所有的光源，且顺序与灯光数组相同。

void DirectLightingIntegrator::Preprocess(const Scene &scene,Sampler &sampler) {
    if (strategy == LightStrategy::UniformSampleAll) {
        // Compute number of samples to use for each light
        for (const auto &light : scene.lights)
            nLightSamples.push_back(sampler.RoundCount(light->nSamples));

        // Request samples for sampling all lights
        for (int i = 0; i < maxDepth; ++i) {
            for (size_t j = 0; j < scene.lights.size(); ++j) {
                sampler.Request2DArray(nLightSamples[j]);
                sampler.Request2DArray(nLightSamples[j]);
            }
        }
    }
}

作为SamplerIntegrator，DirectLightingIntegrator必须实现的主要方式是Li()。其实现的一般形式与WhittedIntegrator::Li()的实现类似：计算与光线相交点处的BSDF，如果表面带有自发光属性，那么则加上自发光的辐射量，以递归的方式跟踪光线跟而实现镜面反射和透射等。

<Compute direct lighting for DirectLightingIntegrator integrator> ≡
if (strategy == LightStrategy::UniformSampleAll)
L += UniformSampleAllLights(isect, scene, arena, sampler,
nLightSamples);
else
L += UniformSampleOneLight(isect, scene, arena, sampler);

为了理解这两种策略的实现，首先回顾一下5.6节介绍的散射方程，它描述了表面上的点p0处ωo方向上的出射辐射度Lo(p,ωo)，由该点处的入射辐射度的全球积分乘以各个方向的BSDF与余弦值得到。,对于DirectLightingIntegrator，我们对子直接从光源发出到物体的辐射量感兴趣，所以我们只需要表现Ld(p,ω)

$L_o(p,\omega_o)=\int_{g^2}f(p,\omega_o,\omega_i)L_d(p,\omega_i)L_d(p,\omega_i)\left|cos\theta_i\right|d\omega_i$

这里可以分解为场景中n个灯光的和：

$\sum^n_{j=1}\int_{g^2}f(p,\omega_o,\omega_i)L_d(p,\omega_i)L_d(j)\left|cos\theta_i\right|d\omega_i$

其中Ld(j)为第N个灯光的发射入射辐射量：

$L_d(p,\omega_i)=\sum_jL_{d(j)}(p,\omega_i)$

一种有效的方法是分别估计式(14.12)中将每一项相加，这是最基本的直接照明策略，也就是UniformSampleAllLights()。UniformSampleAllLights()除了计算直接光照所需的相交点信息和其他参数外，还接受一个handleMedia参数，该参数控制渲染器在计算直接光照时是否应该考虑体积衰减的影响。

Spectrum UniformSampleAllLights(const Interaction &it,
const Scene &scene, MemoryArena &arena, Sampler &sampler,
const std::vector<int> &nLightSamples, bool handleMedia) {
    Spectrum L(0.f);
    for (size_t j = 0; j < scene.lights.size(); ++j) {
        <Accumulate contribution of jth light to L 855>
    }
    return L;
}

这个函数尝试通过Preprocess()，从采样器取得样本数组。

<Accumulate contribution of jth light to L> ≡ 854
const std::shared_ptr<Light> &light = scene.lights[j];
int nSamples = nLightSamples[j];
const Point2f *uLightArray = sampler.Get2DArray(nSamples);
const Point2f *uScatteringArray = sampler.Get2DArray(nSamples);
if (!uLightArray || !uScatteringArray) {
    <Use a single sample for illumination from light 855>
} else {
    <Estimate direct lighting using sample arrays 855>
}

如果所有请求的数组都已被使用，那么代码将通过调用Sampler::Get2D()返回到单个样本估算。

<Use a single sample for illumination from light> ≡ 855
Point2f uLight = sampler.Get2D();
Point2f uScattering = sampler.Get2D();
L += EstimateDirect(it, uScattering, *light, uLight, scene, sampler,
arena, handleMedia);

对于每个灯光样本，EstimateDirect()计算蒙特卡罗估计器对其贡献的值。

<Estimate direct lighting using sample arrays> ≡ 855
Spectrum Ld(0.f);
for (int k = 0; k < nSamples; ++k)
Ld += EstimateDirect(it, uScatteringArray[k], *light, uLightArray[k],
scene, sampler, arena, handleMedia);
L += Ld / nSamples;

在一个具有大量灯光的场景中，不需要总是计算所有被遮挡灯光的直接照明。蒙特卡罗方法提供了一种计算正确的均值结果的方法。以计算两个函数E[f(x)+g(x)]和的期望值为例。如果我们随机计算f(x)或g(x)中的一个，然后将结果乘以2，那么结果的期望值仍然是f(x)+g(x)。这个概念也可以用到任意数目项的和上;这里我们估算直接照明只计算一个随机选择的灯光的结果，并乘以场景灯光数作为结果补偿。

<Integrator Utility Functions> +≡
Spectrum UniformSampleOneLight(const Interaction &it,
const Scene &scene, MemoryArena &arena, Sampler &sampler,
bool handleMedia) {
<Randomly choose a single light to sample, light 856>
Point2f uLight = sampler.Get2D();
Point2f uScattering = sampler.Get2D();
return (Float)nLights *
EstimateDirect(it, uScattering, *light, uLight, scene, sampler,
arena, handleMedia);
}

该方法比UniformSampleOneLight()中使用的均匀方法可能更具创造性。事实上，我们可以任意设置概率，只要我们对结果进行适当的权重调整并且反射点处采样若干灯光概率为非零值。我们在设置这些概率以反映参考点上光对反射光的相对贡献方面做得越好，蒙特卡罗估算器的效率就越高，将方差降低到可接受水平所需的光线就越少。事实上，最终可以取任意数量的样本，只要对它们进行适当的加权。

估算直接光照积分

我们通过估算它的积分，来估算所选的一种特定灯光的直接照明：

$\int_{g^2}f(p,\omega_o,\omega_i)L_d(p,\omega_i)\left|cos\theta_i\right|d\omega_i$

对于灯光项的值，通过选定一个或多个方向ωj，并计算这个估算式子：

$\frac{1}{N}\sum^N_{j=1}\frac{f(p,\omega_o,\omega_i)L_d(p,\omega_j)\left|cos\theta_j\right|}{p(\omega_j)}$

我们将选择方向ωj来进行重要性采样来减少方差。因为BSDF与直接辐射项都较为复杂，很难找到与其乘积匹配的采样分布。在这里，我们将使用基于BSDF的采样分布的一些样本以及一些其他灯光。根据它们各自的特点，这两种抽样方法中的一种可能比另一种有效得多。因此，我们将使用其中更有效的多重重要抽样来减少方差。

如图所示，当BSDF在一组较窄的方向上取较大的值时(该方向比通过对光源进行采样获得的方向集要小得多)，它的采样效率要高得多。另一方面，在相反的情况下，当光源较小且BSDF波瓣的集中度较低时，对光源进行采样会更加有效。

通过多重重要性采样，我们不仅可以同时使用两种抽样方法,而且我们也可以这样做的方式消除了极端的方差的情况下,找到高贡献率的方向。

EstimateDirect()试下了这种方法，并计算一个单独光源样本直接照明的估算值。其handleMedia变量控制是否参与介质衰减影响，其高光参数表示直接光照估计是否考虑了完美高光部分。

Spectrum EstimateDirect(const Interaction &it, const Point2f &uScattering,const Light &light, const Point2f &uLight,const Scene &scene, Sampler &sampler,MemoryArena &arena, bool handleMedia, bool specular) {
    BxDFType bsdfFlags =
        specular ? BSDF_ALL : BxDFType(BSDF_ALL & ~BSDF_SPECULAR);
    Spectrum Ld(0.f);
    // Sample light source with multiple importance sampling
    Vector3f wi;
    Float lightPdf = 0, scatteringPdf = 0;
    VisibilityTester visibility;
    Spectrum Li = light.Sample_Li(it, uLight, &wi, &lightPdf, &visibility);
    VLOG(2) << "EstimateDirect uLight:" << uLight << " -> Li: " << Li << ", wi: "
            << wi << ", pdf: " << lightPdf;
    if (lightPdf > 0 && !Li.IsBlack()) {
        // Compute BSDF or phase function's value for light sample
        Spectrum f;
        if (it.IsSurfaceInteraction()) {
            // Evaluate BSDF for light sampling strategy
            const SurfaceInteraction &isect = (const SurfaceInteraction &)it;
            f = isect.bsdf->f(isect.wo, wi, bsdfFlags) *
                AbsDot(wi, isect.shading.n);
            scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);
            VLOG(2) << "  surf f*dot :" << f << ", scatteringPdf: " << scatteringPdf;
        } else {
            // Evaluate phase function for light sampling strategy
            const MediumInteraction &mi = (const MediumInteraction &)it;
            Float p = mi.phase->p(mi.wo, wi);
            f = Spectrum(p);
            scatteringPdf = p;
            VLOG(2) << "  medium p: " << p;
        }
        if (!f.IsBlack()) {
            // Compute effect of visibility for light source sample
            if (handleMedia) {
                Li *= visibility.Tr(scene, sampler);
                VLOG(2) << "  after Tr, Li: " << Li;
            } else {
              if (!visibility.Unoccluded(scene)) {
                VLOG(2) << "  shadow ray blocked";
                Li = Spectrum(0.f);
              } else
                VLOG(2) << "  shadow ray unoccluded";
            }

            // Add light's contribution to reflected radiance
            if (!Li.IsBlack()) {
                if (IsDeltaLight(light.flags))
                    Ld += f * Li / lightPdf;
                else {
                    Float weight =
                        PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                    Ld += f * Li * weight / lightPdf;
                }
            }
        }
    }

    // Sample BSDF with multiple importance sampling
    if (!IsDeltaLight(light.flags)) {
        Spectrum f;
        bool sampledSpecular = false;
        if (it.IsSurfaceInteraction()) {
            // Sample scattered direction for surface interactions
            BxDFType sampledType;
            const SurfaceInteraction &isect = (const SurfaceInteraction &)it;
            f = isect.bsdf->Sample_f(isect.wo, &wi, uScattering, &scatteringPdf,
                                     bsdfFlags, &sampledType);
            f *= AbsDot(wi, isect.shading.n);
            sampledSpecular = (sampledType & BSDF_SPECULAR) != 0;
        } else {
            // Sample scattered direction for medium interactions
            const MediumInteraction &mi = (const MediumInteraction &)it;
            Float p = mi.phase->Sample_p(mi.wo, &wi, uScattering);
            f = Spectrum(p);
            scatteringPdf = p;
        }
        VLOG(2) << "  BSDF / phase sampling f: " << f << ", scatteringPdf: " <<
            scatteringPdf;
        if (!f.IsBlack() && scatteringPdf > 0) {
            // Account for light contributions along sampled direction _wi_
            Float weight = 1;
            if (!sampledSpecular) {
                lightPdf = light.Pdf_Li(it, wi);
                if (lightPdf == 0) return Ld;
                weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
            }

            // Find intersection and compute transmittance
            SurfaceInteraction lightIsect;
            Ray ray = it.SpawnRay(wi);
            Spectrum Tr(1.f);
            bool foundSurfaceInteraction =
                handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr)
                            : scene.Intersect(ray, &lightIsect);

            // Add light contribution from material sampling
            Spectrum Li(0.f);
            if (foundSurfaceInteraction) {
                if (lightIsect.primitive->GetAreaLight() == &light)
                    Li = lightIsect.Le(-wi);
            } else
                Li = light.Le(ray);
            if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;
        }
    }
    return Ld;
}

<Sample light source with multiple importance sampling>
首先，使用Sample_Li()从灯光采样分布中提取一个样本，它还会返回灯光发出的辐射量以及采样方向的pdf值。

<Compute BSDF or phase function’s value for light sample>
只有当灯光成功采样一个方向并返回非零的辐射量时，EstimateDirect()才会继续计算光线相交计算时的BSDF或者相位函数。

<Compute effect of visibility for light source sample>
如果要考虑参与介质，则从灯光到被找亮点之间的辐射量需要按照光束透射率进行比例计算，这与两点之间参与介质的衰减率有关。

<Add light’s contribution to reflected radiance>
现在就可以计算灯光样本的贡献值。回顾14.2.1节，如果灯光使用delta分布，那么在Sample_Li()和PDF返回的自发光辐射量中都有一个隐含的delta分布项，并且在计算估算的过程中，它们会被抵消。在这种情况下，我们不应该尝试应用多重重要抽样，而应该计算标准估计量。如果这不是一个delta分布光源，那么在ωi方向上采样BSDF的pdf值将会使用BSDF::Pdf()取得。它使用了MIS估算，其权重的计算使用了启发式计算法。

<Sample BSDF with multiple importance sampling>
接下来，使用BSDF采样分布生成一个样本。如果光源的自发光配置包含一个delta分布，那么应该跳过这一步，因为在这种情况下，BSDF采样不可能给出从光源接收到光的方向。否则，可以对BSDF进行采样。

<Account for light contributions along sampled direction wi>
一个重要细节是灯光的pdf与多重重要性采样权重，只有BSDF组件采样ωi方向时其类型为无镜面反射类型时，才会被计算。在镜面反射的情况下，不应该使用MIS，因为灯光采样过程中没有采样镜面反射方向的机会。

<Find intersection and compute transmittance>
给定一个由BSDF或介质的相位函数采样的方向，我们需要知道沿该方向的光线是否与这个特定的光源相交，如果相交，到达着色表面的灯光辐射量是多少？当考虑参与介质时，需要记录光源到交点之间透射率。

<Add light contribution from material sampling>
代码必须考虑到规则形状的面光源以及与之相应的Shape。还有像无限大的面光源这种不含有多边形的光源，它们通过Light::Le()返回样本射线的辐射量。

光照传输公式

光照传输公式(LTE)是描述场景中辐射平衡分布的控制方程。它根据表面自发光属性，给出在表面上的一个点的反射辐射总量、BSDF以及该点的入射照明分布。目前我们只考虑没有参与介质的情况。

一个点的入射光会受到场景内物体的集合形状与散射属性的影响，这会让LTE的评估变得困难。解释这种复杂的渲染算法被称为全局光照算法。区别于光照算法，后者在着色计算中只使用了局部表面属性的信息。

在这一节中，我们将首先推导LTE，并描述一些并叙述一些处理方程的方法，使其更易于求解。

基础理论

LTE的关键原理就是能量守恒。我们必须最终所有的能量变化，因为我们假设照明是一个线性过程，所以系统输入和输出的能量之差也必须等于系统发射和吸收的能量之差。

$\Phi_o-\Phi_i=\Phi_e-\Phi_a$

为了在一个表面上实现能量平衡，那么出色辐射量必须等于自发光加上散射的入射辐射量。自发光辐射量可以用Le表示，散射辐射量可以由散射方程得出，那么Lo为：

$L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{g^2}f(p,\omega_o,\omega_i)L_i(p,\omega_i)\left|cos\theta_i\right|d\omega_i$

因为我们假设没有参与介质，所以沿着场景的光线的辐射是恒定的（碰到任何物体前）。因此，我们可以将p处的入射辐射量与另一点p’处的出射辐射量联系起来。

如果我们定义一个光线投射函数t(p,ω)作为计算光线在ω方向上与表面第一个相交点p的函数。以此我们可以得出p点处入射辐射量与p’点处出射辐射量之间的关系：

$L_i(p,\omega)=L_o(t(p,\omega),-\omega)$

在上述场景处于非封闭的情况下，且光线(p,ω)不与场景中任一物体相交，我们将会定义光线投射函数用于返回特定值 $\Lambda$ ，使得 $L_o(\Lambda,\omega)$ 总是为0。出于简单考虑，这里去掉Lo的下标，所以通过上述关系LTE可以写为：

$L(p,\omega_o)=L_e(p,\omega_o)+\int_{g^2}f(p,\omega_o,\omega_i)L(t(p,\omega_i),-\omega_i)\left|cos\theta_i\right|d\omega_i$

分析LTE的解析方法

LTE的简化方式掩盖了它无法通过解析法来求解。物理BSDF模型的复杂度、任意场景的多边形以及对象间复杂的可见关系交织在一起，因而需要一种有效的数值数值方案对其进行求解。幸运的是，光线追踪算法与蒙特卡洛积分这两个有效的工具可以处理这类复杂的问题，而无需对LTE各组成部分设置限制条件。

这里可以通过简化设置的方式来分析LTE的解析方法，但对通用的渲染方式没有帮助。不过它可以用于积分器的debug。

LTE的表面形式、路径积分、delta分布积分

光看公式都懂了，所以略。

拆分积分表达式

目前已经有许多渲染算法被开发出来，他们在一些特定情况下解LTE十分优秀，但在其他情况却不能很好地工作或根本不能工作。例如Whitted积分器只处理delta分布的BSDF镜面反射，而忽略来自漫反射与光泽的多重散射光线。

因为我们希望能够推导出正确的光线传输算法，并且在不忽略任何贡献与不重复计算其他贡献的情况下计算出所有可能的散射，所以仔细地计算出特解方法在LTE中所占的比例是很重要的。一个解决这个问题的好方法是用不同的方式拆分LTE。例如：我们可以将传输的路径展开：

$L(p_1\rightarrow p_o)=P(\overline{p_1})+P(\overline{p_2})+\sum^\infty_{i=3}P(\overline{p_i})$

其中第一项可以通过计算p1出自发光辐射量获得，第二项使用精确的直接光照解决方案技术来解决，但是其余项将使用更快但不太精确的方法来处理。针对当前渲染场景，如果该附加项对于总体反射辐射度贡献相对较小，则该方法相对合理。

路径追踪

路径追踪是第一个在图形学中使用的通用型无偏差的蒙特卡洛光线传输算法。Kajiya在1986年发表的叙述光线传输方程的论文中介绍了它。路径追踪递增生成
场景中从摄像机到光源的散射事件的路径。从一个方面可以把它看成是whitted方法的扩展，包括了delta分布、非delta分布的BSDFs及光源，而不仅仅是delta项。

虽然直接从基本的光输运方程推导路径跟踪稍微容易一些，但我们将从路径积分形式来解决它，这有助于建立对路径积分方程的理解，并使16.3节介绍的双向路径跟踪更容易理解。

概述

略

路径采样

鉴于这种使用有限项估算无限项的和，我们需要一种估算特定项贡献值的方法 $P(\overline{p_i})$ 。我们需要i+1个顶点来确定路径，其中最后一个顶点pi位于光源上，第一个顶点p0在摄像机上。对于这种 $P(\overline{p_i})$ ,即一个物体表面积的多重积分，其中最自然的事情就是根据场景中物体的表面积生成样本顶点pi，这样与其他样本顶点有着相同的概率。

我们可以定义场景中n个物体的离散概率，如果每个物体的表面积设为Ai，那么第i个物体路径顶点的采样概率为：

$p_i=\frac{A_i}{\sum_jA_j}$

然后，给出了一种对第i个物体上的点进行均匀概率抽样的方法，对物体i上的任意点进行抽样的PDF为1/Ai。因此，采样点的总体概率密度为：

$\frac{A_i}{\sum_jA_j}\frac{1}{A_i}$

所有样本都有着相同的PDF值：

$p_A(p_i)=\frac{1}{\sum_jA_j}$

当给定基于上述采样方法的顶点集合后，可对场景光源上的最后一个顶点pi进行采样，并以同样的方式定义它的pdf。虽然我们可以采用相同技术来对路径顶点到光源采样点进行采样，但这样会产生较大的方差，因为对于所有不在自发光表面的pi点所在路径，其值皆为0。虽然期望值仍为正确的积分，但收敛速度会非常慢。对此一种更好的方法是只对在自发光物体的区域进行采样，并包含更新概率信息。当一条完整路径被确定完毕后，则具备了计算 $P(\overline{p_i})$ 的全部信息。

使用这种通用方法设置采样概率会更具创意。例如：如果我们知道了场景中大部分照明来着若干物体，那么我们就可以为这些物体分配更高的概率来生成路径顶点pi，从而适当地调整采样权重值。

然而这种路径采样方式存在两个相互关联的问题。即产生较大方差或不正确的结果。第一个问题产生的原因是：如果一对相邻的顶点之间不可见，那么他们的贡献值会为0。考虑在复杂的模型中应用这种面积抽样方法:路径中相邻的顶点之间几乎总是有一到两个墙，因此该路径没有贡献，并且在估计中方差很大。

第二个问题产生的原因是：如果被积函数中含有delta分布（点光源或者带有完美反射的BSDF），那么当前采样方案会无法选择路径顶点，使得delte分布包含非0值。同时及时没有delta分布项，随着BSDF的光泽度的增加，所有路径的贡献值都会降低，这将导致bsdf变得更小或出现0值，所以方差将会变大。同样的，未经显式采样的小面积光源也可被视为方差源。

递增式路径构造

针对上述两种问题，一种解决方案是使用递增式路径构造，从摄像机顶点p0开始，在每个顶点采样BSDF以此生成新的方向；下一个顶点pi+1是通过在采样方向上跟踪一条射线并找到最近的交点得到的。我们正在大量的尝试，寻找具有重要的局部贡献的方向，从而有效地寻找具有较大贡献的路径。

因为这个方法是根据立面角采样BSDFs来构造路径的，又因为路径积分LTE是对场景表面积的积分，所以概率密度将根据立体角pω来修正：

$p_A=p_\omega\frac{\left|cos\theta_i\right|}{\mid\mid p_i-p_{i+1}\mid\mid^2}$

因为pi与pi+1互相可见，所以其几何分布项可设为1。

具体实现

路径追踪使用了之前小节介绍的方法，实现了计算路径贡献值 $P(\overline{p_i})$ 之和的估算值。从第一个摄像机光线与场景多边形的交点p1开始，它通过在当前顶点采样BSDF分布，并追踪光线到下一个顶点的方式来递增路径顶点样本。为了找到在特定路径的最后一个顶点，必在场景其中一个光源上的pi，PBRT使用了为方向光积分器开发的基于多重重要性采样的方向光代码。通过使用多重重要性采样 $p_A(p_i)$ 权重，来计算之前描述的估算值，对BSDF采样是寻找光源上的点的更好方法，这结果有着更小的方差。

除了灯光是如何被采样，另一个小区别是在估算路径各项贡献值 $P(\overline{p_i})$ 时，前一个路径顶点pi-1（除开在自发光物体上的点），在构建长度为i的路径时，作为一个起始点是重复利用的。这意味着只需再追踪一条射线即可构建新路径，而不是重头开始构建。那么以这种重复利用路径的方式引入所有具有相关性的 $P(\overline{p_i})$ 进行求和。尽管这将降低当前计算结果的质量，但在实际操作中，这样做减少了需要追踪的光线量，从而提高了渲染效率。

类名为PathIntegrator，继承自SamplerIntegrator。在下面的实现中，L存储了通过路径追踪计算出来的 $P(\overline{p_i})$ 辐射总量，ray存储了追踪下一个路径上顶点的光线。

循环的第一步是通过与场景多边形进行求交测试来寻找下一个路径顶点。
<Intersect ray with scene and store intersection in isect>

如果光线命中一个自发光物体，那物体的自发光属性通常会被忽略，因为在前一个路径顶点的循环迭代已经执行过一次直接光照估算，该估算已经考虑到了它的影响。当光线逃逸到自发光环境中时也是如此。然而，也有两个例外:第一个是与相机光线的初始交点，因为这是直接可见自发光物体的唯一机会。第二是当最后一个路径顶点的采样方向指向镜面反射BSDF组件时：在这种情况下，上一个迭代的直接光照估算不能关联到包含Dirac函数的相关被积函数中。<Possibly add emitted light at intersection>

如果找不到交点，说明光线已经逃离了场景，因此路径迭代终止。当迭代次数超出最大光线反弹次数时，也会终止。
<Terminate path if ray escaped or maxDepth was reached>

当需要计算自发光物体时，在找到交点的情况下或者遇到无线大的面光源发出的辐射量情况下，路径追踪会将权重乘以当前路径顶点的自发光辐射量。
<Add emitted light at path vertex or from the environment>

在估算当前顶点的直接照明前，需要计算当前顶点的散射函数。有一个特殊情况，那就是SurfaceInteraction::bsdf为空指针，此时表明灯光对当前的表面没有影响。PBRT使用这些表面来表示参与介质的边界。由于基础的PathIntegrator忽略了介质，所以它会跳过后续计算。
<Compute scattering functions and skip over medium boundaries>

直接光照计算使用UniformSampleOneLight(),它计算当前路径中的最后一个顶点处的来自直接光照的出射辐射量估算值。
<Sample illumination from lights to find path contribution>

现在很有必要在当前路径的最后一个点采样BSDF，以获得下一条要追踪射线的方向。如上所述，积分器更新路径吞吐量权重，并使用要跟踪的光线初始化光线，以便在for循环的下一次迭代中找到下一个顶点。
<Sample BSDF to get new path direction>

路径追踪将在若干次光线反弹后终止，其终止的概率根据路径吞吐量权重设置。一般来说终止低贡献路径的概率较高，这是值得的，因为它们对最终图像的影响相对较小。
<Possibly terminate the path with Russian roulette>

Spectrum PathIntegrator::Li(const RayDifferential &r, const Scene &scene,Sampler &sampler, MemoryArena &arena,int depth) const {
    ProfilePhase p(Prof::SamplerIntegratorLi);
    Spectrum L(0.f), beta(1.f);
    RayDifferential ray(r);
    bool specularBounce = false;
    int bounces;
    // Added after book publication: etaScale tracks the accumulated effect
    // of radiance scaling due to rays passing through refractive
    // boundaries (see the derivation on p. 527 of the third edition). We
    // track this value in order to remove it from beta when we apply
    // Russian roulette; this is worthwhile, since it lets us sometimes
    // avoid terminating refracted rays that are about to be refracted back
    // out of a medium and thus have their beta value increased.
    Float etaScale = 1;

    for (bounces = 0;; ++bounces) {
        // Find next path vertex and accumulate contribution
        VLOG(2) << "Path tracer bounce " << bounces << ", current L = " << L
                << ", beta = " << beta;

        // Intersect _ray_ with scene and store intersection in _isect_
        SurfaceInteraction isect;
        bool foundIntersection = scene.Intersect(ray, &isect);

        // Possibly add emitted light at intersection
        if (bounces == 0 || specularBounce) {
            // Add emitted light at path vertex or from the environment
            if (foundIntersection) {
                L += beta * isect.Le(-ray.d);
                VLOG(2) << "Added Le -> L = " << L;
            } else {
                for (const auto &light : scene.infiniteLights)
                    L += beta * light->Le(ray);
                VLOG(2) << "Added infinite area lights -> L = " << L;
            }
        }

        // Terminate path if ray escaped or _maxDepth_ was reached
        if (!foundIntersection || bounces >= maxDepth) break;

        // Compute scattering functions and skip over medium boundaries
        isect.ComputeScatteringFunctions(ray, arena, true);
        if (!isect.bsdf) {
            VLOG(2) << "Skipping intersection due to null bsdf";
            ray = isect.SpawnRay(ray.d);
            bounces--;
            continue;
        }

        const Distribution1D *distrib = lightDistribution->Lookup(isect.p);

        // Sample illumination from lights to find path contribution.
        // (But skip this for perfectly specular BSDFs.)
        if (isect.bsdf->NumComponents(BxDFType(BSDF_ALL & ~BSDF_SPECULAR)) >
            0) {
            ++totalPaths;
            Spectrum Ld = beta * UniformSampleOneLight(isect, scene, arena,
                                                       sampler, false, distrib);
            VLOG(2) << "Sampled direct lighting Ld = " << Ld;
            if (Ld.IsBlack()) ++zeroRadiancePaths;
            CHECK_GE(Ld.y(), 0.f);
            L += Ld;
        }

        // Sample BSDF to get new path direction
        Vector3f wo = -ray.d, wi;
        Float pdf;
        BxDFType flags;
        Spectrum f = isect.bsdf->Sample_f(wo, &wi, sampler.Get2D(), &pdf,
                                          BSDF_ALL, &flags);
        VLOG(2) << "Sampled BSDF, f = " << f << ", pdf = " << pdf;
        if (f.IsBlack() || pdf == 0.f) break;
        beta *= f * AbsDot(wi, isect.shading.n) / pdf;
        VLOG(2) << "Updated beta = " << beta;
        CHECK_GE(beta.y(), 0.f);
        DCHECK(!std::isinf(beta.y()));
        specularBounce = (flags & BSDF_SPECULAR) != 0;
        if ((flags & BSDF_SPECULAR) && (flags & BSDF_TRANSMISSION)) {
            Float eta = isect.bsdf->eta;
            // Update the term that tracks radiance scaling for refraction
            // depending on whether the ray is entering or leaving the
            // medium.
            etaScale *= (Dot(wo, isect.n) > 0) ? (eta * eta) : 1 / (eta * eta);
        }
        ray = isect.SpawnRay(wi);

        // Account for subsurface scattering, if applicable
        if (isect.bssrdf && (flags & BSDF_TRANSMISSION)) {
            // Importance sample the BSSRDF
            SurfaceInteraction pi;
            Spectrum S = isect.bssrdf->Sample_S(
                scene, sampler.Get1D(), sampler.Get2D(), arena, &pi, &pdf);
            DCHECK(!std::isinf(beta.y()));
            if (S.IsBlack() || pdf == 0) break;
            beta *= S / pdf;

            // Account for the direct subsurface scattering component
            L += beta * UniformSampleOneLight(pi, scene, arena, sampler, false,
                                              lightDistribution->Lookup(pi.p));

            // Account for the indirect subsurface scattering component
            Spectrum f = pi.bsdf->Sample_f(pi.wo, &wi, sampler.Get2D(), &pdf,
                                           BSDF_ALL, &flags);
            if (f.IsBlack() || pdf == 0) break;
            beta *= f * AbsDot(wi, pi.shading.n) / pdf;
            DCHECK(!std::isinf(beta.y()));
            specularBounce = (flags & BSDF_SPECULAR) != 0;
            ray = pi.SpawnRay(wi);
        }

        // Possibly terminate the path with Russian roulette.
        // Factor out radiance scaling due to refraction in rrBeta.
        Spectrum rrBeta = beta * etaScale;
        if (rrBeta.MaxComponentValue() < rrThreshold && bounces > 3) {
            Float q = std::max((Float).05, 1 - rrBeta.MaxComponentValue());
            if (sampler.Get1D() < q) break;
            beta /= 1 - q;
            DCHECK(!std::isinf(beta.y()));
        }
    }
    ReportValue(pathLength, bounces);
    return L;
}