渲染方程中提到过一个非常重要的概念,双向反射分布函数BRDF,它描述了一个物体在点p处将方向的入射光反射到方向的能力。BRDF是一个简化的模型,只考虑了光的反射,现实中的光照更加复杂,除了反射还有折射、透射等,尽管在程序中没有办法完全模拟出现实世界的光照效果,我们还是希望有一个模型可以尽可能地逼近现实。反射模型(Reflection Models)就是用来完成这个工作的(尽管它还是叫做反射模型),它描述了光照与物体交互的效果,是渲染中不可缺少的一部分。

基本理念

通常来说,反射包括以下四种类型:

  • 漫反射:物体将入射光线均匀地反射到四周。
  • 光泽反射:物体将入射光线向特定的方向散射,从散射方向观察物体可以看到明显的“高光”。
  • 全反射(完美反射):像镜面一样,将入射光线全部反射至反射方向。
  • 逆射:将入射光线沿着入射反向散射射回去,比如天鹅绒材质,会在物体表面形成一层“光泽”。
    物体并不是只进行一种类型的反射,通常是多种反射的混合,但是根据能量守恒定理,反射出的所有光线能量不会大于入射的光线能量。
    光线反射还有一种概念是各向同性和各向异性,对于各向同性的反射,沿着反射点的法线旋转,从不同角度观察到的反射效果是一样的,而对于各向异性的反射,反射的效果会不同。最常见的包含这种反射的物体就是磨砂金属。

反射模型的描述

通常使用BxDF函数来描述反射模型,字母x并没有特殊的含义,是因为几种函数的名称只有第二个字母不一样,BxDF可以分为以下几类:

  • BRDF(双向反射分布函数,Bidirectional Reflectance Distribution Function):用于非透明材质的光照计算。
  • BTDF(双向透射分布函数,Bidirectional Transmission Distribution Function):用于透明材质的光照计算。折射光穿透介质进入另外一种介质时的光照计算模型,只对有透明度的介质适用。
  • BSDF(双向散射分布函数,Bidirectional Scattering Distribution Function):实际上是BRDF和BTDF的综合体。
  • BSSRDF(双向表面散射分布函数,Bidirectional Surface Scattering Reflectance Distribution Function):它常用于模拟透明材质,目前是主流技术。它和BRDF的不同之处在于,BSSRDF可以再现光线透射材质的效果,还可以指定不同的光线入射位置和出射位置。
    各种BxDF函数都是以下形式:

BxDF接口设计

渲染系统中BxDF的接口描述了BxDF的通用函数,具体的设计如下:

class BxDF
    : 
{
  public:
       BxDFFlags Flags() const;
       using TaggedPointer::TaggedPointer;
       
       std::string ToString() const;
       SampledSpectrum f(Vector3f wo, Vector3f wi, TransportMode mode) const;
       pstd::optional<BSDFSample>
       Sample_f(Vector3f wo, Float uc, Point2f u,
                TransportMode mode = TransportMode::Radiance,
                BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;
       Float PDF(Vector3f wo, Vector3f wi, TransportMode mode,
                 BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;
       SampledSpectrum rho(Vector3f wo, pstd::span<const Float> uc,
                           pstd::span<const Point2f> u2) const;
       SampledSpectrum rho(pstd::span<const Point2f> u1,
                                  pstd::span<const Float> uc2,
                                       pstd::span<const Point2f> u2) const;
       void Regularize();
};

下面是每个函数的含义。

BxDFFlags Flags()

返回此物体BxDF标志,它由物体的材质决定,BxDFFlags的定义为:

enum BxDFFlags {
    Unset = 0,
    Reflection = 1 << 0,
    Transmission = 1 << 1,
    Diffuse = 1 << 2,
    Glossy = 1 << 3,
    Specular = 1 << 4,
       DiffuseReflection = Diffuse | Reflection,
       DiffuseTransmission = Diffuse | Transmission,
       GlossyReflection = Glossy | Reflection,
       GlossyTransmission = Glossy | Transmission,
       SpecularReflection = Specular | Reflection,
       SpecularTransmission = Specular | Transmission,
       All = Diffuse | Glossy | Specular | Reflection | Transmission
};

Reflection和Transmission决定了光线是进行反射还是透射,下面的Diffuse,Glossy和Specular决定了反射(透射)的类型,这里将逆射并入了Glossy类型。

f函数

这是BxDF的核心函数,它接受出射光线方向wo,入射光线方向wi以及光线传输模式mode,返回BxDF函数的值。
按照BxDF函数的定义,f函数应该还需要反射发生的具体的位置点p,但是在程序设计中,必然是先找到了光线与物体的点才能找到该点的材质,进而找到使用BxDF,所以这里就不需要点p作为入参。
它的返回值SampledSpectrum部分定义如下:

static constexpr int NSpectrumSamples = 4;
 
class SampledSpectrum {
  public:
    ... 
  private:
    pstd::array<Float, NSpectrumSamples> values;
};

它是一个数组,每个元素包含了4个float值,数组长度代表了采样光谱的数量,后续可以用它们对光谱进行采样,其中每一个元素可以看作一个颜色。
mode参数用来指示出射方向是朝向相机还是朝向光源,或者是入射方向的相反方向,后续可以用来处理不对称的散射,它的定义如下:

enum class TransportMode { Radiance, Importance };

Sample_f函数

Sample_f函数用来提供在已知光线出射方向时的入射方向及其pdf,这个说法看上去有点奇怪,其实在光线追踪算法中,我们是按照光线传输的逆方向发射光线的,光线在人眼(相机)处产生,根据光路的可逆性,最终计算的结果与正常的光照相同。这样做产生一个后果,我们是先知道了出射方向的光线(它是相机发射光线的反方向),然后再去求入射光线的方向。
Sample_f函数的作用很明显,BxDF的函数值是在渲染方程积分项之中的,对它的入射光线方向采样必然要使用到函数分布的pdf。在现实中这种方式也是非常合理的,当物体进行全反射时(比如镜面),除了出射方向的反射方向,其他方向的入射光线是不提供任何共享的,它们的pdf应该等于0。
Sample_f函数定义如下:

pstd::optional<BSDFSample>
Sample_f(Vector3f wo, Float uc, Point2f u,
         TransportMode mode = TransportMode::Radiance,
         BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;

wo为出射方向,uc和u是随机值,比如可以使用uc来决定进行透射还是反射,使用u来选定入射方向(极坐标系),mode上面已经介绍过,sampleFlags强制指定进行透射还是进行反射,或者二者都有。

enum class BxDFReflTransFlags {
    Unset = 0,
    Reflection = 1 << 0,
    Transmission = 1 << 1,
    All = Reflection | Transmission
};

函数的返回值为pstd::optional<BSDFSample>,这表示采样是有可能失败的,比如在不透明的物体上进行透射。
BSDFSample包含了此次采样的结果,它的定义如下,

struct BSDFSample {
    // 构造函数省略
    SampledSpectrum f;    // 此次采样的结果,BxDF的函数值
    Vector3f wi;          // 入射方向
    Float pdf = 0;        // 此次采样立体角的pdf值
    BxDFFlags flags;      // 样本特征,上面有介绍相关定义
    Float eta = 1;        // IOR(折射率)
    bool pdfIsProportional = false;  // 某些随机光路的BxDF会使用到,大多数情况下为false
};

同样,有一些工具函数来快速判断采样的样本特征,

bool IsReflection() const { return pbrt::IsReflective(flags); }
bool IsTransmission() const { return pbrt::IsTransmissive(flags); }
...

PDF函数

PDF函数返回BxDF指定出射和入射方向的pdf值,

Float PDF(Vector3f wo, Vector3f wi, TransportMode mode,
          BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const;

mode参数和sampleFlags参数的含义与Sample_f相同。

rho

这个函数是半球面反射的相关函数,使用蒙特卡洛积分计算物体表面的某个点的辐射率,它有两种类型的重载。
第一种函数提供出射方向wo和采样点数对应的样本,计算spectrum,

SampledSpectrum BxDF::rho(Vector3f wo, pstd::span<const Float> uc,
                          pstd::span<const Point2f> u2) const {
    SampledSpectrum r(0.);
    for (size_t i = 0; i < uc.size(); ++i) {
           pstd::optional<BSDFSample> bs = Sample_f(wo, uc[i], u2[i]);
           if (bs)
               r += bs->f * AbsCosTheta(bs->wi) / bs->pdf;
 
    }
    return r / uc.size();
}

uc和u2是样本值,uc是一维的样本,u2是二维的样本。
其中使用的函数AbsCosTheta:

Float AbsCosTheta(Vector3f w) { return std::abs(w.z); }

因为向量是定义在局部空间中的,该点的法线作为坐标系的z轴,而是计算该向量与法线的夹角,所以向量z坐标乘以向量长度即为,由于wi是归一化的单位向量,所以直接取z坐标的值即可。这里取得是z坐标的绝对值,说明向量位于法线负半球时,与正半球得到的结果相同。
使用公式来表示:

第二种函数并不直接提供wo,而是随机在法线的半球面选择一个wo方向,用于计算所有方向的入射光都相同时半球面的反射率,

SampledSpectrum BxDF::rho(pstd::span<const Point2f> u1,
        pstd::span<const Float> uc, pstd::span<const Point2f> u2) const {
    SampledSpectrum r(0.f);
    for (size_t i = 0; i < uc.size(); ++i) {
           Vector3f wo = SampleUniformHemisphere(u1[i]);
           if (wo.z == 0)
               continue;
           Float pdfo = UniformHemispherePDF();
           pstd::optional<BSDFSample> bs = Sample_f(wo, uc[i], u2[i]);
           if (bs)
               r += bs->f * AbsCosTheta(bs->wi) * AbsCosTheta(wo) / (pdfo * bs->pdf);
 
    }
    return r / (Pi * uc.size());
}

可以看到wo是从SampleUniformHemisphere函数采样获得的。SampleUniformHemisphere用于在半球面上均匀采样,它的实现为:

Vector3f SampleUniformHemisphere(Point2f u) {
    Float z = u[0];
    Float r = SafeSqrt(1 - Sqr(z));
    Float phi = 2 * Pi * u[1];
    return {r * std::cos(phi), r * std::sin(phi), z};
}

这里传入一个二维随机数u,生成半球面上的三维向量(x, y, z),具体的实现原理先忽略,就是根据蒙特卡洛积分和采样中的重要性采样方法推导出来的。
UniformHemispherePDF是半球面上均匀采样的PDF,它是一个常数

Float UniformHemispherePDF() { return Inv2Pi; }

这个函数使用公式来表示:

表示BxDF在入射光在所有方向都相同时,整个半球表面的反射率。

BSDF接口设计

BSDF接口是BxDF接口的简单封装,包含了从渲染空间(可能是世界空间也可能是相机空间)转换的局部空间的变换方法,因为通常BSDF的计算在局部空间完成。

class BSDF {
  public:
       ...
  private:
       BxDF bxdf;
       Frame shadingFrame;
 
};

BSDF并不继承BxDF接口,而是包含一个BxDF成员变量,shadingFrame保存当前渲染帧的信息。
下面是BSDF接口中包含几个重要的方法。

构造函数

BSDF() = default;
BSDF(Normal3f ns, Vector3f dpdus, BxDF bxdf)
    : bxdf(bxdf),
      shadingFrame(Frame::FromXZ(Normalize(dpdus), Vector3f(ns))) {}

接受着色点的法线和bxdf变量,以及dpdu(目前只关心法线即可)。

坐标转换

Vector3f RenderToLocal(Vector3f v) const { return shadingFrame.ToLocal(v); }
Vector3f LocalToRender(Vector3f v) const {
    return shadingFrame.FromLocal(v);
}

用来进行渲染空间与局部空间的相互转换。

f函数

SampledSpectrum f(Vector3f woRender, Vector3f wiRender,
                  TransportMode mode = TransportMode::Radiance) const {
    Vector3f wi = RenderToLocal(wiRender), wo = RenderToLocal(woRender);
    if (wo.z == 0) return {};
    return bxdf.f(wo, wi, mode);
}

计算指定出射和入射方向的BSDF函数值,将向量变换到局部空间后,调用bxdf的接口进行计算。
它还要另一种模版形式,在指定了具体的bxdf类型后直接调用该类型的f函数,

template <typename BxDF>
SampledSpectrum f(Vector3f woRender, Vector3f wiRender,
                  TransportMode mode = TransportMode::Radiance) const {
    Vector3f wi = RenderToLocal(wiRender), wo = RenderToLocal(woRender);
    if (wo.z == 0) return {};
    const BxDF *specificBxDF = bxdf.CastOrNullptr<BxDF>();
    return specificBxDF->f(wo, wi, mode);
}

Sample_f函数

pstd::optional<BSDFSample> Sample_f(
        Vector3f woRender, Float u, Point2f u2,
        TransportMode mode = TransportMode::Radiance,
        BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const {
    Vector3f wo = RenderToLocal(woRender);
    if (wo.z == 0 ||!(bxdf.Flags() & sampleFlags)) return {};
       pstd::optional<BSDFSample> bs = bxdf.Sample_f(wo, u, u2, mode, sampleFlags);
       if (!bs || !bs->f || bs->pdf == 0 || bs->wi.z == 0)
           return {};
       bs->wi = LocalToRender(bs->wi);
       return bs;
}

同样也是调用bxdf的Sample_f函数,这里增加了坐标系的转换,以及对采样结果的检查,如果采样失败或者pdf为0,亦或是采样的wi与法线垂直(bswi.z == 0),那么不会返回采样结果。

PDF函数

Float PDF(Vector3f woRender, Vector3f wiRender,
          TransportMode mode = TransportMode::Radiance,
          BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const {
    Vector3f wo = RenderToLocal(woRender), wi = RenderToLocal(wiRender);
    if (wo.z == 0) return 0;
    return bxdf.PDF(wo, wi, mode, sampleFlags);
}

检测wo方向后调用bxdf的PDF函数。

漫反射(Diffuse Reflection)

漫反射是最简单的一种反射模型,它假设物体表面接受到光照后,会向四面八方均匀地反射。通常这种模型被称为Lambertian反射模型,它的BSDF函数是一个常数。即

假设物体的反射率为R(通常这个反射率就是物体的基础颜色),根据能量守恒定理:

可以解得

所以

DiffuseBxDF的具体设计:

class DiffuseBxDF {
  public:
    ...
  private:
    SampledSpectrum R;
};

R为物体的反射率(反照率),通常是物体的基础颜色。

f函数

SampledSpectrum f(Vector3f wo, Vector3f wi, TransportMode mode) const {
    if (!SameHemisphere(wo, wi))
        return SampledSpectrum(0.f);
    return R * InvPi;
}

如果wo和wi在同一半球,则返回,否则返回0。

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

Sample_f函数

pstd::optional<BSDFSample> Sample_f(
        Vector3f wo, Float uc, Point2f u, TransportMode mode,
        BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const {
    if (!(sampleFlags & BxDFReflTransFlags::Reflection))
        return {};
    Vector3f wi = SampleCosineHemisphere(u);
    // 保证wo和wi在同一半球
    if (wo.z < 0) 
       wi.z *= -1;
    Float pdf = CosineHemispherePDF(AbsCosTheta(wi));
    return BSDFSample(R * InvPi, wi, pdf, BxDFFlags::DiffuseReflection);
}

通过SampleCosineHemisphere函数采样出方向wi,这里的SampleCosineHemisphere函数是对半球上的项进行采样,采样算法与上面的SampleUniformHemisphere相同,同样是重要性采样方法。然后使用CosineHemispherePDF计算在半球面上的pdf。
按理来说,根据渲染方程中的描述,这里的采样只涉及到BxDF项,项是不需要考虑的,但是实践证明,将项加入BxDF项一起考虑会得到更佳的渲染效果,所以这里采样的函数其实是。后面的各类BSDF和BRDF模型也会遵循此原则。

PDF函数

Float PDF(Vector3f wo, Vector3f wi, TransportMode mode,
          BxDFReflTransFlags sampleFlags = BxDFReflTransFlags::All) const {
    if (!(sampleFlags & BxDFReflTransFlags::Reflection) ||
        !SameHemisphere(wo, wi))
        return 0;
    return CosineHemispherePDF(AbsCosTheta(wi));
}

根据上面的原则,这里也是返回CosineHemispherePDF(AbsCosTheta(wi))。

镜面反射和透射的原理

IOR

IOR(Index of Refraction)意味折射率,它反应了光线进入物体时,光速的衰减情况,比如一个材质的IOR为2,那就说明光线进入这个材质后,传播速度会变为光速的一半。通常情况下材质的IOR在1.0到2.5之间。
IOR通常使用来表示。
下面是一些常见介质的IOR。

介质Index of refraction(
真空1.0
1.31
水(20摄氏度)1.333
玻璃1.5-1.6
钻石2.42

镜面反射定律

假设入射光的角度为(,),当它发生镜面反射时,其反射方向为入射方向绕法线旋转180度,使用公式来表示:

如果使用立体角来表示,反射方程推导过程如下。
假设入射方向为,反射方向与入射方向、法线必然位于同一平面,在这个平面上,将分解:

则有:

根据镜面反射定律,有

最终反射方向

使用代码实现:

Vector3f Reflect(Vector3f wo, Vector3f n) {
    return -wo + 2 * Dot(wo, n) * n;
}

折射定律

当入射光以角度为(,)照射物体并发生折射时,折射光线方向与两边介质的IOR相关,具体为

这就是折射定律。
为了方便,通常使用相对折射率来表示折射定律,相对折射率定义为

因此有

在多数情况下,光线追踪算法的光线是从真空射入物体的,所以,此时,计算更加方便。
同样,如果使用立体角来表示,如下图

则有

最终折射方向

代码实现:

bool Refract(Vector3f wi, Normal3f n, Float eta, Float *etap, Vector3f *wt) {
    Float cosTheta_i = Dot(n, wi);
    // cos theta小于0表示光线从物体向外发生折射
    if (cosTheta_i < 0) {
        eta = 1 / eta;
        cosTheta_i = -cosTheta_i;
        n = -n;
    }
 
    Float sin2Theta_i = std::max<Float>(0, 1 - Sqr(cosTheta_i));
    Float sin2Theta_t = sin2Theta_i / Sqr(eta);
    // sin2Theta_t >= 1表示发生全反射
    if (sin2Theta_t >= 1)
        return false;
 
    Float cosTheta_t = SafeSqrt(1 - sin2Theta_t);
 
    *wt = -wi / eta + (cosTheta_i / eta - cosTheta_t) * Vector3f(n);
    // 根据需要提供eta
    if (etap)
        *etap = eta;
 
    return true;
}

菲涅尔方程

上面的反射定律和折射定律描述了光线发生镜面反射和折射时出射光的方向,现在需要知道出射光的能量大小。
当光线照射到一个物体上时,通常会同时发生反射和折射,菲涅尔方程描述了光线的反射率,即光线有多少能量被反射。
菲涅尔方程与入射光和折射光的振幅有关,其中平行方向和垂直方向的计算方式也不相同。而振幅与光的折射率和入射角相关,具体的方程为:

按照之前的表示方式,令。则上面的方程可以表示为:

最终的菲涅尔反射率由两个方向的振幅比值结合起来:

这个比例就是光线发生反射时反射光线所占的能量比例。
代码实现:

Float FrDielectric(Float cosTheta_i, Float eta) {
    cosTheta_i = Clamp(cosTheta_i, -1, 1);
    // cosTheta_i 小于0时反转朝向
    if (cosTheta_i < 0) {
        eta = 1 / eta;
        cosTheta_i = -cosTheta_i;
    }
    // 计算 cosTheta_t
    Float sin2Theta_i = 1 - Sqr(cosTheta_i);
    Float sin2Theta_t = sin2Theta_i / Sqr(eta);
    if (sin2Theta_t >= 1)
        return 1.f;
    Float cosTheta_t = SafeSqrt(1 - sin2Theta_t);
 
    Float r_parl = (eta * cosTheta_i - cosTheta_t) /
                   (eta * cosTheta_i + cosTheta_t);
    Float r_perp = (cosTheta_i - eta * cosTheta_t) /
                   (cosTheta_i + eta * cosTheta_t);
    return (Sqr(r_parl) + Sqr(r_perp)) / 2;
}

导体的菲涅尔方程

导体发生反射和折射行为时,它的IOR是一个复数,通常使用实部描述光速的下降,而虚部描述光线在导体内的衰减,光线在导体内的衰减很快,甚至可以影响到反射分量。使用复数并不是IOR真的具有复数的性质,而是将其作为一个数学工具。
光线的振幅在物体中的衰减可以表示为,其中i为虚数单位,为频率,为折射率IOR,z为距离。假设现在物体的折射率为依然是一个实数,根据复数三角函数的性质,

那么函数E可以表示为
如果IOR为一个复数,其中依然表示折射率,k表示衰减程度,那么函数E可以表示为

不论使用哪种表示方式,虚数单位i都是没有使用到的。
使用这种方式来表示菲涅尔方程,具体的实现为

Float FrComplex(Float cosTheta_i, pstd::complex<Float> eta) {
    using Complex = pstd::complex<Float>;
    cosTheta_i = Clamp(cosTheta_i, 0, 1);
    Float sin2Theta_i = 1 - Sqr(cosTheta_i);
    // 由于eta是一个复数,所以计算出来的sin2Theta_t也是一个复数
    Complex sin2Theta_t = sin2Theta_i / Sqr(eta);
    Complex cosTheta_t = pstd::sqrt(1 - sin2Theta_t);
 
    Complex r_parl = (eta * cosTheta_i - cosTheta_t) /
                     (eta * cosTheta_i + cosTheta_t);
    Complex r_perp = (cosTheta_i - eta * cosTheta_t) /
                     (cosTheta_i + eta * cosTheta_t);
    // pstd::norm计算复数的模
    return (pstd::norm(r_parl) + pstd::norm(r_perp)) / 2;
}

为了方便使用,对其进行一个简单的封装,输入eta和k并对其Fr进行评估。

SampledSpectrum FrComplex(Float cosTheta_i, SampledSpectrum eta,
                          SampledSpectrum k) {
    SampledSpectrum result;
    for (int i = 0; i < NSpectrumSamples; ++i)
        result[i] = FrComplex(cosTheta_i,
                              pstd::complex<Float>(eta[i], k[i]));
    return result;
}

粗糙度和微平面理论

现实中物体发生反射时,完美的镜面反射是不存在的,因为没有完美的光滑镜面,物体越粗糙,它的镜面反射就会越弱。
从微观上说,物体的表面是由许多凹凸不平的微平面组成的,光线照射到物体的表面时,这些微平面同时发生反射,物体的表面越光滑(此时微平面之间的法线差异越小,即微平面都朝向同一个方向),整体的反射效果就越强,而物体的表面越粗糙(微平面朝向不同的方向),整体的反射效果就越弱。这就是微平面理论的基础。
按照上面这个理论,其实可以在模型上使用许多细小的三角形来模拟真实的粗糙表面,然后在渲染时对每个微平面来计算反射效果,这样也可以表示出不同的表面粗糙度,但是这对于建模和光线追踪算法来说都是十分困难和低效的,所以通常情况下,会使用统计近似的方法来对微平面理论进行建模。

微平面的法线分布函数