PBRT阅读：第十一章 纹理 第11.5-11.6节

http://www.opengpu.org/forum.php?mod=viewthread&tid=6034

11.5 立体纹理和过程纹理

一旦我们将用于2D纹理函数的(s,t)纹理坐标视为可以由任意函数生成的量（而不仅仅被视为表面的参数坐标），我们就很自然地将纹理函数进一步一般化，将之定义为三维域上的函数（常称为立体纹理）。立体纹理的一个方便之处在于，所有物体自然而然地有一个三维纹理映射---即利用物体空间位置进行纹理映射。对于那些没有二维参数化形式的物体（如三角形网格，隐式方程表达的曲面），还有参数形式产生纹理变形的参数化物体（如球面上两极上的变形），三维纹理是很有优势的。第11.2.5节为此做了一些准备工作，我们定义了通用的TextureMapping3D接口来计算3D纹理坐标，还具体实现了类IdentityMapping3D。

然而立体纹理会产生一个新问题：纹理的表示问题。三维的图像贴图要占据相当大的存储空间，而且跟直接可以从照片或美术作品获取的二维纹理贴图相比，获取三维图像要难得多。因此，在开发立体纹理技术的过程中，过程纹理的技术就应运而生了：该技术使用短小的程序在表面上的任何位置生成纹理值。一个过程纹理的例子是一个过程式的正弦波。假定我们用一个正弦波来做凹凸贴图（例如，模拟水面的波纹），如果在一个网格上的点上预先计算好函数值并存到纹理图中，就很没有效率，而且有潜在的误差。实际上，在需要的时候直接在表面上求sin()函数值更有意义。

举例来说，如果我们找到描述木块内纹理的颜色的三维函数，就可以生成任何形似木雕的复杂物体的图像。许多年来，过程纹理技术已经可以过程性地描述越来越复杂的表面，过程纹理在实际应用中有了相当的进步。

过程纹理有几个有趣的事实。首先，由于不需要高分辨率的纹理贴图，也就可以了减少内存需求。还有，过程式的着色可允许无限量的细节；当观察者靠近物体时，纹理函数可以在着色点上求值，从而很自然地生成可被观察到的足够多的细节。相反地，当观察者离物体太近时，图像贴图就显得模糊不清。然而，过程纹理在对纹理外观的细微之处的控制方面，要比图像贴图要困难得多。

另一个过程纹理的困难之处是反走样。过程纹理的求值常常是很费时的，而且过程纹理并不具备图像贴图（在反走样过程中）十分有效的采样点集。因为我们期望在采样之前就去掉纹理函数中的高频信息，所以我们就必须在每个步骤都要留意高频内容，从而避免高频的介入。虽然这听起来很不容易，但有不少很好用的技术来解决这个问题。

11.5.1 UV纹理

我们的第一个过程纹理将表面(u,v)坐标转换为Spectrum的前两个分量。这个纹理在对新形体类的程序纠错时很有用。


<UVTexture Declarations> =
class UVTexture : public Texture<Spectrum> {
    public:
        <UVTexture Public Methods>
    private:
        TextureMapping2D *mapping;
    };

<UVTexture Public Methods> =
    Spectrum Evaluate(const Differential Geometry &dg) const {
        float s, t , dsdx, dtdx, dsdy, dtdy;
        mapping->Map(dg, &s, &t, &dsdx, &dtdx, &dsdy, &dtdy);
        float cs[COLOR_SAMPLES];
        memset(cs, 0, COLOR_SAMPLES * sizeof(float));
        cs[0] = s - Floor2Int(s);
        cs[l] = t - Floor2Int(t);
        return Spectrum(cs);
    }

11.5.2 棋盘纹理

棋盘纹理是一种很标准的过程纹理。我们用(s,t)纹理坐标将参数空间划分成正方形区域，并且这些区域的着色是用交替的样式进行的。这里我们没有用两个固定的颜色作为棋盘格的样式，而是允许用户传入两个纹理来对这些区域交替着色。我们可以传入两个ConstantTexture来得到传统的黑白棋盘。


<CheckerboardTexture Declarations> =
    template <class T> class Checkerboard2D : public Texture<T>{
    public:
        <Checkerboard2D Public Methods>
    private:
        <Checkerboard2D Private Data>
    };

为了方便，我们设格子函数(check function)在(s,t)空间中的频率为1： 即在每个方向上的格子的宽度为一个单位长。当然，我们可以用带有适当的(s,t)比例因子的TextureMapping2D类可以对之进行调整。

<Checkerboard2D Public Methods> =
    Checkerboard2D(TextureMapping2D *m, Reference<Texture<T> > c1,
            Reference<Texture<T> > c2, const string &aa) {
        mapping = m;
        tex1 = c1;
        tex2 = c2;
        <Select antialiasing method for Checkerboard2D>
    }   
   
Checkerboard2D Private Data =
    Reference<Texture<T> > tex1, tex2;
    TextureMapping2D *mapping;

使用棋盘纹理可以很好地验证对各种过程纹理反走样方法的权衡效果。这里我们实现了三种反走样方式，并用一个字符串来表示。

<Select antialiasing method for Checkerboard2D> =
    if (aa == "none") aaMethod = NONE;
    else if (aa == "super-sample") aaMethod = SUPERSAMPLE;
    else if (aa == "closedform") aaMethod = CLOSEDFORM;
    else {
        Warning("Anti-aliasing mode \"%s\" not understood "
            "by Checkerboard2D, defaulting"
            "to \"supersample\"", aa.c_str());
        aaMethod = SUPERSAMPLE;
    }

<Checkerboard2D Private Types> =
    enum { NONE, SUPERSAMPLE, CLOSEDFORM } aaMethod;

求值例程做通常的纹理坐标和微分计算，然后使用相应的程序来计算反走样的棋盘值。

<Checkerboard2D Public Methods> +=
    T Evaluate(const Differential Geometry &dg) const {
        float s, t , dsdx, dtdx, dsdy, dtdy;
        mapping->Map(dg, &s, &t, &dsdx, &dtdx, &dsdy, &dtdy);
        if (aaMethod == CLOSEDFORM) {
            <Compute closed-form box-filtered Checkerboard2D value>
        }
        else if (aaMethod == SUPERSAMPLE) {
            <Supersample Checkerboard2D>
        }

        <Point sample Checkerboard2D>
    }

最简单的情形是忽略反走样，只是对棋盘纹理做点采样。对于这个情形，在用TextureMapping2D得到(s,t)纹理坐标之后，就可以计算出(s,t)所对应的整数的棋盘坐标，将两个棋盘x,y坐标值加在一起，通过检查结果的奇偶位来确定是用哪一个纹理来求值。

<Point sample Checkerboard2D> =
    if ((Floor2Int(s) + Floor2Int(t) % 2 == 0)
        return texl->Evaluate(dg);
    return tex2->Evaluate(dg);

由于对棋盘纹理做点采样的效果走样很严重，我们必须做点努力进行适当的反走样。最简单的情况是整个滤波区域落在一个棋盘格之内。这时我们只需确定格子的纹理类型并对之求值。只要格子的纹理本身是反走样的，其结果也应该是反走样的。

<Compute closed-form box-filtered Checkerboard2D value> =
    <Evaluate single check if filter is entirely inside one of them>
    <Apply box filter to checkerboard region>

我们可以先计算滤波区域的包围盒，如果包围盒的范围是在棋盘格之内，则整个滤波区域就落在棋盘格之内。本节将使用由偏导数(∂s/∂x, ∂s/∂y,等等)计算出来的和轴对齐的包围盒来做滤波，而不是使用前面所介绍的EWA滤波器中的椭圆区域。这样就简化了实现，当然也会增加了滤波值的模糊度。在下面的代码中，变量ds和dt分别是在两个方向上的滤波器宽度的一半，所以整个滤波范围是从(s-ds, t - dt)到(s+ds, t+dt)。

<Evaluate single check if filter is entirely inside one of them> =
    float ds = max(fabsf(dsdx), fabsf(dsdy));
    float dt = max(fabsf(dtdx), fabsf(dtdy));
    float s0 = s - ds, s1 = s + ds;
    float t0 = t - dt, t1 = t + dt;
    if (Floor2Int(s0) == Floor2Int(s1) && Floor2Int(t0) == Floor2Int(t1)){
        <Point sample CheckerboardZD>
    }

否则的话，我们就计算一个浮点数比率来指明滤波区域覆盖两种格子类型的情况。这等价于计算2D阶梯函数(在tex1中取0，在tex2中取1)在整个滤波区域中的平均值。如图中的(a)表示这个阶梯函数，(b)是对之进行从0到x的积分：
   

    棋盘函数c(x)和相应的积分函数定义如下：

   
有了平均值之后，就可以根据这个比率将两个纹理进行混合。 为了计算阶梯函数在两个方向上的平均值，我们分别在两个方向上做棋盘函数的积分：

<Apply box filter to checkerboard region> =
    #define BUMPINT(x) \
        (Floor2Int((x)/2) + \
        2.f * max((x/2)-Floor2Int(x/2) - .5f, 0.f))
    float sint = (BUMPINT(s1) - BUMPINT(s0)) / (2.f * ds);
    float tint = (BUMPINT(t1) - BUMPINT(t0)) / (2.f * dt);
    float area2 = sint + tint - 2.f * sint * tint;
    if (ds > 1.f || dt > 1.f)
        area2 = .5f;
    return (1.f - area2) * tex1->Evaluate(dg) +
            area2 * tex2->Evaluate(dg);

最后我们要实现的反走样方法是超采样。该方法在纹理空间中在点(s,t)附近取一组随机分层(random stratified)的位置对棋盘函数求值，这些位置经过抖动，使得它们可以大致可以覆盖滤波区域。超采样的一个有利之处是前面所介绍的解析式盒滤波器实际上并不完全正确。特别是，它假定通过确定滤波器覆盖两种格子类型的比例并将两种纹理做相应的混合来确定总的反走样结果。这个假定所带来的问题是，每个子纹理都各自求值并进行反走样，好像它们在整个滤波区域都是可见的。

下图是显示这个错误假定而人为制作的例子：如果两个子纹理本身是棋盘，按照如图的结构做计算，则整个棋盘的纹理变成了灰色。


当然这是人为制造出的极端例子，对于只是常量的子纹理而言，盒滤波器方法还是可以正确工作的。然而，超采样方法没有这个问题。

超采样在s,t方向上分别使用固定数目的采样（N_SAMPLES）。对于每个采样位置，以原DifferentialGeometry对象为初始值计算一个新的DifferentialGeometry对象，并对棋盘纹理进行求值。

<Supersample Checkerboard2D> =
    #define SQRT_SAMPLES 4
    #define N_SAMPLES (SQRT_SAMPLES * SQRT_SAMPLES)
    float samples[2*N_SAMPLES];
    StratifiedSample2D(samples, SQRT_SAMPLES, SQRT_SAMPLES);
    T value = 0.;
    float filterSum = 0 . ;
    for (int i = 0; i < N_SAMPLES; ++i) {
        <Compute new differential geometry for supersample location>
        <Compute (s, t) for supersample and evaluate subtexture>
    }
    return value / filterSum;

对于滤波区域中的每个采样点，该采样的DifferentialGeometry可以由(s,t)处的偏导数来近似地求得，其参数坐标也可以通过采样的偏移量计算出来。

<Compute new differential geometry for supersample location> =
    float dx = samples[2*i] - 0.5f;
    float dy = samples[2*i+l] - 0.5f;
    DifferentialGeometry dgs = dg;
    dgs.p += dx * dgs.dpdx + dy;
    dgs.u += dx * dgs.dudx + dy;
    dgs.v += dx * dgs.dvdx;
    dgs.dudx /= N_SAMPLES;
    dgs.dudy /= N_SAMPLES;
    dgs.dvdx /= N_SAMPLES;
    dgs.dvdy /= N_SAMPLES;

最后，这些采样值要通过一个高斯滤波器做加权平均，从相应子纹理计算出的值被累加到value中。

<Compute (s, t) for supersample and evaluate subtexture> =
    float ss, t s , dsdxs, dtdxs, dsdys, dtdys;
    mapping->Map(dgs, &ss, &ts, &dsdxs, &dtdxs, &dsdys, &dtdys);
    float wt = expf( -2.f * (dx*dx + dy * dy));
    filterSum += wt;
    if ((Floor2Int(ss) + Floor2Int(ts)) % 2 == 0)
        value += wt * tex1->Evaluate(dgs);
    else
        value += wt * tex2->Evaluate(dgs);

11.5.3 立体棋盘纹理

前面介绍的Checkerboard2D类在参数空间中将物体包上一个棋盘纹理样式。我们也可以定义一个基于三维纹理坐标的立体棋盘纹理，使得物体看上去是由三维棋盘格子雕刻出来似的。象二维棋盘纹理一样，这个实现根据一个查找函数来确定要使用的纹理函数。需要注意的是这两个纹理本身不需要是立体纹理，Checkerbord3D根据点的三维位置来选择纹理。如图是一个使用了三维棋盘纹理的例子：
   

<CheckerboardTexture Declarations> +=
    template <class T> class Checkerboard3D : public Texture<T> {
    public:
        <Checkerboard3D Public Methods>
    private:
        <Checkerboard3D Private Data>
    };


<CheckerboardsD Public Methods> =
    Checkerboard3D(TextureMapping3D *m,
            Reference<Texture<T> > c1,
            Reference<Texture<T> > c2)
        mapping = m;
        texl = c1;
        tex2 = c2;
    }

<Checkerboard3D Private Data> =
    Reference<Texture<T> > tex1 , tex2;
    TextureMapping3D *mapping;

如果忽略反走样，确定一个点是否在一个3D格子里的区域的基本计算如下：

    ((Floor2Int(P.x) + Floor2Int(P.y) + Floor2Int(P.z)) % 2 == 0).

对于反走样而言，这里的实现使用了跟2D棋盘中的超采样方式。这里略去了代码，因为除了判定点是否在格子里的方法不同以外，其它都是相同的。

<Checkerboard3D Public Methods> +=
    T Evaluate(const Differential Geometry &dg) const {
        <Supersample Checkerboard3D>
    }

11.6 噪声函数

为了写出可以描述复杂表面外观的立体纹理程序，我们介绍一些很有用的对立体纹理的生成过程进行控制的机制。
设想有一个由许多木板条组成的地板，每个木板条的颜色都有些不同。或者考虑一下一个被风吹拂的湖面，我们想在这个湖面上呈现具有类似波幅的水浪，但并不希望所有的水浪是均匀(homogeneous)的（这是当所有波是一组正弦波的叠加时所发生的情况）。在纹理中对这类变化的模拟有助于提高最终结果的真实性。

开发这类纹理的一个困难之处在于渲染器在一组不规则分布的点上进行纹理求值，每个求值过程都是彼此独立的。就这样，过程纹理隐式地(implicitly)定义了一个复杂的样式，对该样式上的所有值的查询只限于对这些点上的值的计算。与此相反的是显式(explict)样式的描述，例如PostScript语言，它将一页上的图形用一系列的绘图命令来描述。隐式的方式有一个困难之处：即纹理不能在每个求值点上仅仅靠调用RandomFloat()来得到随机性，因为这样做的效果使得每个点有与它的邻点全然不同的随机值，对于所生成的纹理样式而言没有任何连贯性。

关于这个问题（即对图形学中的过程纹理引入可控的随机性）有一个优雅的解决方案，即使用所谓的噪声函数(noise function)。一般而言，图形学中的噪声函数将n维空间(n=1,2,3)平滑地映射到[-1,1]而没有明显的重复。一个实际应用中的噪声函数的最关键的性质是它有已知的最大频率且是有带宽局限的。这样我们就可以控制噪声函数加入到纹理的频率内容，使得高于Nyquist极限所允许的频率无法进入纹理。

许多噪声函数是基于三维空间中的整数格(integer lattice)开发出来的。首先，空间中的每个整数位置(x,y,z)对应一个值，然后我们对这些格点值(lattice values)进行插值来计算特定点上的噪声值。这个思想即可一般化为更高的维数或限制到更低的维数d，这样一个格单元的点格数就是2^d。这个方法的一个简单的例子是值噪声(value noise)，即每个格点上对应一个-1到1之间的伪随机数，而实际的噪声值是通过三线性插值或更复杂的样条插值计算出来的，这样就可以避免格单元之间的导数不连续性，从而得到更平滑的结果。

对于这样的噪声函数，必须能够高效地计算出给定整数(x,y,z)格点上的参数值。因为存贮所有(x,y,z)点上的值是不可能的，所以需要一些更紧凑的表示方法。其中一个选择是使用一个散列表，用坐标做键值来从一个大小固定的表格中查找预先计算好的参数值。

11.6.1 Perlin噪声函数

在pbrt中，我们实现一个由Ken Perlin发明的噪声函数，即Perlin噪声函数。它在所有(x,y,z)的整数格点上的值都是0。其变化是由每个格点上的梯度向量控制的，这些向量用来指导一个平滑函数在格点之间的插值。这个噪声函数有许多前面所介绍的优良特性，它在计算上是高效的，且易于实现。
   

<Texture Method Definitions> +=
    COREDLL float Noise(float x, float y, float z) {
        <Compute noise cell coordinates and offsets>
        <Compute gradient weights>
        <Compute trilinear interpolation of weights>
    }

为了方便，我们有下列一个使用一个点的噪声函数：

<Texture Method Definitions> +=
    COREDLL float Noise(const Point &P) {
        return Noise(P.x, P.y, P.z);
    }

首先我们计算给定点所在的格单元的整数坐标，以及该点距离单元左下角的小数偏移量：

<Compute noise cell coordinates and offsets> =
    int ix = Floor2Int(x);
    int iy = Floor2Int(y);
    int iz = Floor2Int(z);
    float dx = x - ix, dy = y - iy, dz = z - iz;

然后我们计算8个单元角所对应的权值。每个整数格点有一个梯度向量与之对应。每个角上的梯度向量对单元内任何一点的影响是通过计算梯度向量和单元角(corner)到给定点的向量之间的点积来得到的，这个计算由函数Grad()来处理。如图，实线是梯度向量，虚线是角到给定点的向量。



<Compute gradient weights> =
    ix &= (NOISE_PERM_SIZE-1);
    iy &= (NOISE_PERM_SIZE-1);
    iz &= (NOISE_PERM_SIZE-1);
    float w000 = Grad(ix, iy, iz, dx, dy, dz);
    float w100 = Grad(ix+1, iy, iz, dx-1, dy, dz);
    float w010 = Grad(ix, iy+1, iz, dx, dy-1, dz);
    float w110 = Grad(ix+1, iy+1, iz, dx-1, dy-1, dz);
    float w001 = Grad(ix, iy, iz+1, dx, dy, dz-1);
    float w101 = Grad(ix+1, iy, iz+1, dx-1, dy, dz-1);
    float w011 = Grad(ix, iy+1, iz+1, dx, dy-1, dz-1);
    float w111 = Grad(ix+1, iy+1, iz+1, dx-1, dy-1, dz-1);

我们通过对一个预计算的整数值表格(NoisePerm)进行索引来查找某个特定的格点上的梯度向量。格点值的前4位确定了使用16个梯度向量中的哪一个。在预处理阶段，我们将表格的NOISE_PERM_SIZE个表项填上0到NOISE_PERM_SIZE-1，然后再随机地将其排列。然后复制这些值，做成一个2*NOISE_PERM_SIZE大小的数组，即原表格内容再被重复一次：这样做的目的是提高下列代码的效率：

给定了特定的(ix,iy,iz)格点以后，用下列方法来查值：

    NoisePerm[NoisePerrn[NoisePerm[ix]+iy]+iz];

用这个方法来查值（而不是用诸如NoisePerm[ix+iy+iz]的方法）可以是结果更具不规则性，例如，当ix和iy互换后，该方法不会返回相同的值。另外，由于表格被复制后增大了一倍，就省去了下面代码的取模计算：

    (NoisePerm[ix]+iy) % NOISE_PERM_SIZE

通过查表确定了梯度向量号之后，下一步是计算相应的点积。然而，我们无需梯度向量的直接表示。所有的梯度向量用-1,0,1这三个值为分量，所以点积可以被简化为向量的分量相加。最后的代码实现如下：

<Texture Method Definitions> +=
    inline float Grad(int x, int y, int z, float dx, float dy, float dz) {
        int h = NoisePerm[NoisePerm[NoisePerm[x]+y]+z];
        h &= 15;
        float u = h<8 || h==12 || h==13 ? dx : dy;
        float v = h<4 || h==12 || h==13 ? dy : dz;
        return ((h&1) ? -u : u) + ((h&2) ? -v : v);
    }

<Perlin Noise Data> =
    #define NOISE_PERM_SIZE 256
    static int NoisePerm[2 * NOISE_PERM_SIZE] = {
        151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96,
        53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142,
        <Remainder of the noise permutation table>
    };

有了梯度向量的8个贡献值之后，下一步是在点上做三线性插值。我们并不是用dx,dy,dz直接插值，而是将其传入一个平滑函数NoiseWeight()。这样可以确保噪声函数在格单元之间有一阶和二阶导数连续。

<Texture Method Definitions> +≡
    inline float NoiseWeight(float t) {
        float t3 = t*t*t;
        float t4 = t3*t;
        return 6.f*t4*t - 15.f*t4 + 10.f*t3;
    }

<Compute trilinear interpolation of weights> =
    float wx = NoiseWeight(dx), wy = NoiseWeight(dy), wz = NoiseWeight(dz);
    float x00 = Lerp(wx, w000, w100);
    float x10 = Lerp(wx, w010, w110);
    float x01 = Lerp(wx, w001, w101);
    float x11 = Lerp(wx, w011, w111);
    float y0 = Lerp(wy, x00, x10);
    float y1 = Lerp(wy, x01, x11);
    return Lerp(wz, y0, y1);

11.6.2 随机波尔卡圆点(Random Polka Dots)

噪声函数的一个基本用途是用来生成波尔卡圆点纹理，这种纹理将(s,t)纹理空间划分为长方形的单元，每个单元内有圆点的几率为50%，圆点被随机地放置到它们所在的单元内。DotsTexture使用一个2D映射函数，还有两个纹理，其中一个纹理用于圆点之外的区域，另一个用于圆点之内的区域。
   

<DotsTexture Declarations> =
    template <class T> class DotsTexture : public Texture<T> {
    public:
        <DotsTexture Public Methods>
    private:
        <DotsTexture Private Data>
    };

<DotsTexture Public Methods> =
    DotsTexture(TextureMapping2D *m, Reference<Texture<T> > c1,
            Reference<Texture<T> > c2) {
        mapping = m;
        outsideDot = c1;
        insideDot = c2;
    }
<DotsTexture Private Data> =
    Reference<Texture<T> > outsideDot, insideDot;
    TextureMapping2D *mapping;

求值函数首先使用(s,t)纹理坐标计算sCell和tCell的整数值，即所在的单元坐标。这里没有考虑波尔卡圆点纹理的反走样，留做练习。

<DotsTexture Public Methods> +=
    T Evaluate(const DifferentialGeometry &dg) const {
        <Compute cell indices for dots>
        <Return insideDot result if point is inside dot>
        return outsideDot->Evaluate(dg);
    }
<Compute cell indices for dots> =
    float s, t , dsdx, dtdx, dsdy, dtdy;
    mapping->Map(dg, &s, &t, &dsdx, &dtdx, &dsdy, &dtdy);
    int sCell = Floor2Int(s + .5f) , tCell = Floor2Int(t + .5f);

有了单元坐标以后，就需要判定该单元是否有波尔卡圆点。显然，我们需要这个计算有一致性，使得单元内所有的点都能产生同样的结果。然而，我们却希望结果是不规则的（例如，我们不想要隔一个单元出现一个圆点的情况）。噪声函数可以解决这个问题：对于所有单元内的点用同一个位置(sCell+0.5, tCell+0.5)来求噪声函数的值，这样的结果即有一致性又有不规则性。如果该值大于零，就在单元内放入一个圆点。

如果单元内有圆点，我们再利用噪声函数对圆点的中心进行随机地移位。在对该点进行噪声函数求值时，先对之偏移任意的常量，这样该点所使用的噪声值来自不同的噪声单元，这样做是为了去掉前面判定单元是否有圆点的噪声求值所引起的相关性。

<Return insideDot result if point is inside dot> =
    if (Noise(sCell+.5f, tCell+.5f) > 0) {
        float radius = .35f;
        float maxShift = 0.5f - radius;
        float sCenter = sCell + maxShift *
            Noise(sCell + 1.5f, tCell + 2.8f);
        float tCenter = tCell + maxShift *
            Noise(sCell + 4.5f, tCell + 9.8f);
        float ds = s - sCenter, dt = t - tCenter;
            if (ds*ds + dt*dt < radius*radius)
        return insideDot->Evaluate(dg);
    }

10.6.3 噪声函数的特性和光谱合成

噪声函数是有带宽局限的，这一事实意味着我们可以通过对定义域做比例变换来调整它的频率内容。例如，如果Noise(p)有已知的频率内容，那么Noise(2*p)的频率内容就可以提高一倍。这正如sin(x)和sin(2x)之间的频率内容的关系一样。我们可以利用这项技术来创建其变化率符合我们需要的噪声函数。

在过程纹理的许多应用中，受多个比例因子控制的变化特别有用。例如，我们可以在基本噪声函数上加上更精致的变化。其中一个有效的方法是通过光谱合成(spectral synthesis)来计算纹理样式，即一个复合函数fs(x)是由另一个函数f(x)通过一组贡献值的加权和来定义的：
        
其中wi是权值，si是参数值的比例因子。如果基函数f(x)有良好定义的频率内容(例如，它是一个正弦或余弦函数，或者是一个噪声函数），那么每一项f(six)也有良好的频率内容。因为叠加和中的每一项有一个权值wi，则其结果就是不同频率的贡献值的叠加，并且具有不同权值的频率范围。

在一般情况下，我们用几何级数来选择比例因子s_i，使得s_i = 2 s_i-1，并且权值满足w_i = w_i-1 / 2。其效果就是当我们加入更高的频率变化时，它对整个f_s(x)的总体形状有较小的影响。每个添加项被称为噪声函数的一个倍频程(octave)，因为它的频率内容是前一项的两倍。当我们对Perlin噪声函数使用这个规定是，就得到所谓的方形布朗运动fBm(fractional Brownian motion)。
分形布朗运动是非常有用的过程纹理的创建工具，因为它比普通的噪声函数有更复杂的变化，容易计算，而且有良好定义的频率内容。工具函数FBm()实现了分形布朗运动。

除了要进行函数求值的点以及点处的偏导数之外，该函数还用到一个参数omega，它的范围是0到1，它通过控制高频贡献值的衰减程度来影响纹理样式的平滑度，还有maxOctaves，它是求和计算是所用的最大倍频程的数目。

<Texture Method Definitions> +=
    float FBm(const Point &P, const Vector &dpdx, const Vector &dpdy,
            float omega, int maxOctaves) {
        <Compute number of octaves for antialiased FBm>
        <Compute sum of octaves of noise for FBm>
        return sum;
    }

FBm函数中的反走样基于“截取”(clamping)技术，它的思想是：我们通过对一组分量进行累加而得到结果值，其中每个分量具有已知的频率内容，当分量的加入导致出现高于Nyquist极限的频率时，就要停止加入分量。因为Noise()函数的平均值为零，我们只需计算出其中没有过高频率的倍频程的个数，而不需要对更高的倍频程进行函数求值。

Noise()（fs(x)的第一项也是如此）的最大频率内容为ω = 1，则每个后继项有加倍的频率内容。因此，我们希望找到适当的前n项，使得：

        2ⁿs = 2ω = 2

其中s为噪声空间中的采样速率。这个条件保证了在给定的采样速率下不会有无法表示的频率。

这样，我们有：

        2ⁿ-1 = 1/s
        n - 1 = log (1/s)
        n = 1 - 1/2 log(s²)

采样速率的平方s²可以通过计算偏导数∂p/∂x和 ∂p/∂y的长度中较大的一个来得到：

<Compute number of octaves for antialiased FB> =
    float s2 = max(dpdx.LengthSquared(), dpdy.LengthSquared());
    float foctaves = min((float)maxOctaves, 1.f - .5f * Log2(s2));
    int octaves = Floor2Int(foctaves);

最后octaves（不超过Nyquist极限的）的整数部分所对应的倍频程被累加在一起，最后一个倍频程也根据foctaves的小数部分也加入进来。这样就保证后续的倍频程能够渐进地加入，而不是突然地出现，在过渡区出现显著的人为缺陷。这里增加倍频程之间的频率的倍数是1.99（而不是2），这样会减小噪声函数在格点为0所产生的影响。

<Compute sum of octaves of noise for FBm> =
    float sum = 0., lambda = 1., o = 1.;
    for (int i = 0; i < octaves; ++i) {
        sum += o * Noise(lambda * P);
        lambda *= 1.99f;
        o *= omega;
    }
    float partialOctave = foctaves - octaves;
    sum += o * SmoothStep(.3f, .7f, partialOctave) * Noise(lambda * P);

SmoothStep()函数的参数包括一个最大值和一个最小值，以及一个要做平滑插值的点。如果点在最小值之下，则返回0，如果在最大值之上，则返回1。否则，用一个三次Hermit样条函数进行插值计算：

<Texture Inline Functions> =
    inline float SmoothStep(float min, float max, float value) {
        float v = Clamp((value - min) / (max - min), 0.f, 1.f);
        return v * v * (-2.f * v + 3.f);
    }

跟FBm()密切相关的是Turbulence()函数。它也是对噪声函数组成的若干项进行叠加，但是取每一项的绝对值：
   
取绝对值的运算造成了合成函数的一阶导数不连续。这样turbulence函数有无限的频率内容。然而，这个函数的视觉特征很适用于过程纹理，如图：

turbulence()试图使用FBm()的方法来反走样。但由于一阶导数不连续，其结果并非完全成功。但是这个反走样过程至少去掉了某些最坏的情况；当然最后的求助方法是增大像素采样速率。在实际应用中，这个函数在过程纹理中的走样并非无可救药，特别是跟几何边界和阴影边界的无限高频比较而言。

<Texture Method Definitions> +=
    float Turbulence(const Point &P, const Vector &dpdx, const Vector &dpdy,
            float omega, int maxOctaves) {
        <Compute number of octaves for antialiased FBm>
        <Compute sum of octaves of noise for turbulence>
        return sum;
    }


<Compute sum of octaves of noise for turbulence> =
    float sum = 0., lambda = 1., o = 1.;
    for (int i = 0; i < octaves; ++i) {
        sum += o * fabsf(Noise(lambda * P));
        lambda *= 1.99f;
        o *= omega;
    }

    float partialOctave = foctaves - octaves;
    sum += o * SmoothStep(.3f, .7f, partialOctave) *
        fabsf(Noise(lambda * P));

11.6.4 凹凸纹理和皱褶纹理

fBM()和turbulence()函数对于凹凸贴图而言非常有用，因为它们可以在纹理中加入随机变化。FBmTexture是一个浮点数纹理，它用FBm()计算偏移量，而WrinkedTexture是用turbulence()计算。
        
<FBmTexture Declarations> =
    template <typename T> class FBmTexture : public Texture<T> {
    public:
        <FBmTexture Public Methods>
    private:
        <FBmTexture Private Data>
    };

<FBmTexture Public Methods> =
    FBmTexture(int oct, float roughness, TextureMapping3D *map)
        : omega(roughness), octaves(oct), mapping(map) { }
   
<FBmTexture Private Data> =
    float omega;
    int octaves;
    TextureMapping3D *mapping;

<FBmTexture Public Methods> +=
    T Evaluate(const DifferentialGeometry &dg) const {
        Vector dpdx, dpdy;
        Point P = mapping->Map(dg, &dpdx, &dpdy);
        return FBm(P, dpdx, dpdy, omega, octaves);
    }

WrinkedTexture的实现跟FBmTexture几乎一样，只是用turbulence()调用代替FBm（），这里省略了。


16.5.5 风中的波浪

使用fBm可以产生令人信服的波浪效果。这个纹理基于两个观察。首先，在一个被风吹拂的湖面，有些区域相对平缓，有些却很汹涌，这是因为风在不同区域的自然变化而产生的。再者，表面上单个波浪可以用基于fBm并用风的强度进行缩比的波浪样式来描述。

<WindyTexture Declarations> =
    template <typename T> class WindyTexture : public Texture<T> {
    public:
        <WindyTexture Public Methods>
    private:
        <WindyTexture Private Data>
    };

<WindyTexture Public Methods> =
    WindyTexture(TextureMapping3D *map) : mapping(map) { }

<WindyTexture Private Data>
    TextureMapping3D *mapping;


这里的求值函数调用了两次FBm()函数。 第一个调用将点P进行10倍的缩比(即乘以0.1)；这样FBm()返回相对较低的频率变化。这个值用来确定风的局部强度。第二个调用确定波浪的跟风无关的波幅。这两个值的乘积就是波浪的实际偏移量。

<WindyTexture Public Methods> +=
    T Evaluate(const DifferentialGeometry &dg) const {
        Vector dpdx, dpdy;
        Point P = mapping->Map(dg, &dpdx, &dpdy);
        float windStrength = FBm(.1f * P, .1f * dpdx, .1f * dpdy, .5f, 3);
        float waveHeight = FBm(P, dpdx, dpdy, .5f, 6);
        return fabsf(windStrength) * waveHeight;
    }

16.5.6 大理石

噪声函数的另一个典型用途是在运用另一个纹理或查找表之前对纹理坐标进行扰动。例如，我们可以用这样的方法模拟大理石：将大理石看作一系列的分层，然后用噪声函数扰动用来查找分层的坐标值。本节的MarbleTexture实现了这个方法。下图显示了这个纹理背后的思想：左边的图直接使用了球面上y的坐标对大理石的分层进行查找；右图则先对y坐标值进行fBM扰动，从而产生了纹理上的变化。
   
<MarbleTexture Declarations> =
    class MarbleTexture : public Texture<Spectrum> {
    public:
        <MarbleTexture Public Methods>
    private:
        <MarbleTexture Private Data>
    };

该纹理的构造器使用了FBm()函数所用的那些参数。另外还有参数variation，用来指定扰动的幅度。


<MarbleTexture Public Methods> =
    MarbleTexture(int oct, float roughness, float sc, float var,
        TextureMapping3D *map)
    : octaves(oct), omega(roughness), scale(sc), variation(var),
    mapping(map) { }

<MarbleTexture Private Data> =
    int octaves;
    float omega, scale, variation;
    TextureMapping3D *mapping;

对给定点的y值扰动后，就得到该点在大理石分层中的偏移量，再用正弦函数将之映射到[0,1]之内。然后，<Evaluate marbel spline at t>使用t值对一个三次样条函数（由跟真大理石颜色相近的一系列颜色定义）求值。

<MarbleTexture Public Methods> +=
    Spectrum Evaluate(const DifferentialGeometry &dg) const {
        Vector dpdx, dpdy;
        Point P = mapping->Map(dg, &dpdx, &dpdy);
        P *= scale;
        float marble = P.y + variation *
            FBm(P, scale * dpdx, scale * dpdy, omega, octaves);
        float t = .5f + .5f * sinf(marble);
        <Evaluate marble spline at t>
    }