PBRT阅读：第五章 颜色和辐射度学 第5.1节

第五章 颜色和辐射度学

为了生成图像，我们需要精确地描述光和对光的采样方法，这就要求我们必须具备一些辐射度学的背景知识。辐射度学的研究对象是在某一环境中电磁辐射的传播。其中我们感兴趣的是用于渲染过程中的可见光，即波长在370nm到730nm之间的电磁波。低端的波长（大约400nm）对应于蓝光，中端的波长（大约550nm）对应于绿光，高端的波长（大约650nm）对应于红光。

在本章，我们将介绍四个描述电磁辐射的量：辐射通量(flux)，辐射强度(intensity)，辐射照度(irradiance)和辐射亮度(radiance)。这些量都可以由它们的光谱功率分布(spectral power distribution, SPD)来定义。SPD是关于波长的分布函数，它描述了每个波长上的光的总量值。第5.1节定义的Spectrum类就是来表示SPD的。

5.1 光谱表示

现实世界中的物体的光谱功率分布是很复杂的。下图分别是一个日光灯发出的光和一个柠檬皮的反射光的光谱功率分布图。利用SPD进行计算的渲染器需要一个简洁、有效和精确的方法来表示这些分布函数。在具体实现中，当然要做质量上的权衡。


关于这类问题的研究通常起始于寻找能够表示SPD的良好的基函数。基函数的意义在于把无限空间中的SPD函数映射到关于几个实数系数ci的低维空间中。例如，一个简单的基函数是常量函数B(λ) = 1。任何SPD函数都可以由一个值等于其平均值的单个系数c来表示，这样它的近似表示就是cB(λ) =c。当然，这是一个质量不高的近似，因为绝大多数的SPD比这个单一的基函数所能表示的SPD要复杂得多。

5.1.1 Spectrum类

在pbrt中，Specturm类及其内建的操作符实现了所有关于SPD的计算。Spectrum类隐藏了光谱表示的细节，如果要改变系统的这类细节，只需改变Spectrum类的实现，而无需改变其它的代码。Specturm类的声明在core/color.h，并定义在core/color.cpp。

然而，我们没有使用插件的方法，不同的Spectrum实现并不以插件的方式跟系统对接；如果改变了具体的实现，整个系统就要重新编译。这个设计的好处是很多Spectrum类中的函数是用inline写的，速度比用虚函数的方式要快得多。把经常使用的短函数写成内联(inline)函数会使系统的性能有实质性的提高。

这里所讨论的Spectrum的标准实现基于最直接了当的表示方法：它存储在固定波长上采样的一组采样值。这个表示方法对于变化相对平缓SPD还是很好的，也可以做为评判其它表示方法的基准。如果该方法采用大量的采样值（例如按1nm的间隔采样），就可以作为一种参考实现（当然这要花费大量的计算资源）。这种方法的基本操作（例如SPD的加运算和乘运算）还是很高效的，其复杂性跟采样值的个数呈线性关系。

这种表示的缺点是不能适应那些含有尖峰的SPD（比如日光灯的SPD）。如果采样值错过了尖峰，就会产生非常不准确的图像。
Spectrum类存放采样值的个数是由宏定义COLOR_SAMPLES在编译时确定的。

<Global Constants> +=
    #define COLOR_SAMPLES 3
<Spectrum Declarations> =
    class COREDLL    Spectrum {
    public:
        <Spectrum Public Methods>
        <Spectrum Public Data>
    private:
        < Spectrum Private Data>
    };
< Spectrum Private Data> =
    float c[COLOR_SAMPLES];

我们提供了两个构造器函数，第一个用同一个值初始化所有的采样值，另一个用采样值数组进行初始化：
<Spectrum Public Methods> =
    Spectrum(float v = 0.f) {
        for (int i = 0; i < COLOR_SAMPLES; ++i)
            c[ i ] = v;
    }

<Spectrum Public Methods> +=
    Spectrum(float cs[COLOR_SAMPLES]) {
        for (int i = 0; i < COLOR_SAMPLES; ++i)
            c[ i ]= cs[ i ];
    }

Spectrum类还支持算术操作，实现也很简单。第一个操作是光谱分布的叠加。很容易证明两个SPD的和的每个采样值等于两个SPD相应采样值的和。

<Spectrum Public Methods> +=
    Spectrum &operator += (const Spectrum &s2) {
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            c[ i ] += s2.c[ i ];
        return *this;
    }

<Spectrum Public Methods> +=
    Spectrum operator + (const Spectrum &s2)  const {
        Spectrum ret = *this;
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            ret.c[ i ]  += s2.c[ i ];
        return ret;
    }

类似地我们可以实现Spectrum的减、乘、除操作，它们的实现跟加法操作类似，从略。

由于pbrt有一段性能要求很高的代码，即要用一个Spectrum对象的加权值来修正另一个Spectrum对象，所以我们要引入能够更有效地实现这个操作的函数AddWeighted()。（更底层的原因是，对于象s = s + w*s2这类的表达式，一些编译器无法摆脱不必要的暂时变量）。

<Spectrum Public Methods> +=
    void AddWeighted(float w, const Spectrum &s) {
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            c[ i ]  += w * s.c[ i ];
    }

我们还提供了一个测试相等的函数：

<Spectrum Public Methods> +=
    bool operator == (const Spectrum &sp) const {
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            if(c [ i ] != sp.c[ i ]) return false;
        return true;
    }

有时我们需要知道一个Spectrum对象是否代表一个全零的SPD。例如，如果一个表面完全不反射光，光传输例程就可以不生成反射光线，从而避免了额外的计算开销。

<Spectrum Public Methods> +=
    bool Black () const {
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            if(c[ i ]  != 0.) return false;
        return true;
    }

Spectrum还包括几个更少见的几个操作函数，包括取平方根函数和幂函数。它们会在第九章中的Fresnel类和Lafortune表面反射模型中用到。

<Spectrum Public Methods> +=
    Spectrum Sqrt () const {
        Spectrum ret;
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            ret.c [ i ] = sqrtf(c[ i ]);
        return ret;
    }

<Spectrum Public Methods> +=
    Spectrum Pow (const Spectrum &e) const {
        Spectrum ret;
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            ret.c [ i ]= c [ i ] > 0 ? powf(c[ i ], e.c[ i ]) : 0.f;
        return ret;
    }

在第12,17章有关参与介质的计算中，我们还会用到取负值操作和指数操作。

<Spectrum Public Methods> +=
    Spectrum operator - () const ；
    friend COREDLL Spectrum Exp(const Spectrum &s);

在图像处理管线中，我们有时会把Spectrum值限制在给定范围值内：

<Spectrum Public Methods> +=
    Spectrum Clamp (float low = 0.f, float hight = INFINITY) const {
        Spectrum ret;
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            ret.c [ i ] = ::Clamp(c[ i ], low, hight);
        return ret;
    }
最后，我们提高一个出错辅助函数，用来检查SPD的采样值中是否有NaN(Not a Number)值。当出现被零除的情况时，就会有这种现象发生。

<Spectrum Public Methods> +=
    bool IsNaN () const {
        for (int i = 0; i < COLOR_SAMPLES; ++i) 
            if(isnan(c[ i ])) return false;
        return true;
    }

5.1.2 XYZ颜色

人的视觉系统的一个显著的特性是视网膜上的视锥细胞只对红色、绿色和蓝色最敏感，这样我们就可以用三个浮点数来表示颜色。关于颜色感知的三刺激理论说明了所有可见的SPD都可以精确地用三个值X_λ, Y_λ，Z_λ来表示。给定一个SPD(λ)，这些值可以由它和光谱匹配曲线函数X(λ), Y(λ), Z(λ)的乘积的积分而得到。
              
这三个曲线是由国际照明委员会(CIE)标准组织通过大量的试验而得到的（如图）。我们可以认为这些匹配曲线类似于视网膜上的对三类颜色（红，绿，蓝）产生感知的视锥细胞的反应曲线。值得注意的是，分布情况很不相同的SPD可能有非常相接近的X_λ, Y_λ，Z_λ值。对于人类观察者而言，这样的SPD的视觉效果是一样的。这样的光谱分布被称为条件等色(metamers)。
 

这就给我们带来一些关于光谱功率分布表示上的小麻烦。大多数颜色空间试图描述可见色，因而只用三个系数，充分地利用了感知颜色的三刺激理论。虽然XYZ可以很好地表示对于人类观察者所看到的SPD，但对于光谱计算而言，这并不是一组很好的基函数。例如，XYZ值可以很好地表示柠檬皮上的颜色和日光灯上的颜色，但是它们的XYZ乘积所产生的XYZ颜色跟求得更精确的SPD乘积后再计算XYZ值的XYZ颜色很不一致。

了解了这一点限制后，我们要在Spectrum类中加一个求XYZ值的函数。这个函数会被ImageFilm类用到，它可以把最后的用Spectrum表示的图像转换成适合于显示的RGB颜色。在这个过程中，第一步就是先把Spectrum转换成跟显示方法无关(display-indepentent)的XYZ值。

对于点采样的光谱表示，用来表示X_λ, Y_λ，Z_λ的SPD和XYZ匹配函数的卷积实际上就是对点采样值的加权求和。

<Spetrum Public Methods> +=
    void XYZ(float xyz[3]) const {
        xyz[0] = xyz[1] = xyz[2] = 0.;
        for(int i = 0; i < COLOR_SAMPLES; ++i) {
            xyz[0] += XWeight * c[ i ];
            xyz[1] += YWeight * c[ i ];
            xyz[2] += ZWeight * c[ i ];
        }
    }

因此，我们还要确定SPD采样点所对应的波长。如果绕过这个问题来实现Spectrum类也是可能的。其中一个权益之计就是使用CRT彩色显像管中的标准红绿蓝三原色的光谱，虽然对于高质量的光谱计算而言这是不够的。高清晰电视标准已经有了一组关于这些RGB光谱的标准定义。把这些RGB值转换成XYZ值的权值如下：

<Spectrum Method Definitions> =
    float Spectrum::XWeight[COLOR_SAMPLES] = {
        0.412453f, 0.357580f, 0.180423f
    };
    float Spectrum::YWeight[COLOR_SAMPLES] = {
        0.212671f, 0.715160f, 0.072169f
    };
    float Spectrum::XWeight[COLOR_SAMPLES] = {
        0.019334f, 0.119193f, 0.950227f
    };

我们还提供另一个有用的创建Spectrum对象的函数，它把XYZ值转换成内部表示。对于点采样的Spectrum表示而言，这仍然是一个加权和：

<Spectrum Method Definitions> +=
    Spectrum FromXYZ(float x, float y, float z) {
        float c[3];
        c[0] = 3.240479f * x + -1.537150f * y + -0.498535f * z;
        c[1] = -0.969256f * x + 1.875991f * y + 0.041556f * z;
        c[2] = 0.055648f * x + -0.204043f * y + 1.057311f * z;
        return Spectrum(c);
    }

为了方便地计算出其它光谱基函数所对应的XWeight,YWeight,ZWeight值，我们提供了从360nm到830nm按每1nm间隔采样的标准X(λ), Y(λ), Z(λ)反应曲线上的值。

<Spectrum Public Data> =
    static const int CIEStart = 360;
    static const int CIEEnd = 360;
    static const int nCIE = CIEend - CIEstart + 1;
    static const float CIE_X[nCIE];
    static const float CIE_Y[nCIE];
    static const float CIE_Z[nCIE];

XYZ颜色中的y值跟光亮度(luminance)密切相关，它描述了一种颜色的明亮程度（照度会在第8.4.1节有更详细的讨论）。为了方便起见，我们提供了计算y的函数：

<Spectrum Public Methods> +=
    float y() const {
        float v = 0.;
        for (int i = 0; i < COLOR_SAMPLES; ++i)
            v += YWeight * c[ i ];
        return v;
    }

我们还可以用光亮度对Spectrum实例从暗到亮进行排序：

<Spectrum Public Methods> +=
    bool operator<(const Spectrum &s2) const {
        return y() <s2.y();
    }