PBRT阅读：第七章 采样和重构 第7.5.2-7.6节

http://www.opengpu.org/forum.php?mod=viewthread&tid=5148&fromuid=10107

http://www.opengpu.org/bbs/forum.php?mod=viewthread&tid=5161&fromuid=10107

7.5.2 使用最佳候选者模式

BestCandidateSampler使用预先生成好的采样表，相当简单。如果每个像素平均有pixelSamples个采样，那么单个采样表的拷贝可以在x或y方向上覆盖SQRT_SAMPLE_TABLE_SIZE/sqrt(pixelSamples)个像素。跟StratifiedSampler和LDSampler一样，该采样器从裁剪窗口的左上角开始，从左到右然后从上到下地扫描图像。在这里，它生成采样表之内的所有采样之后，再前进到采样表所要覆盖的下一个图像区域。下图表示了这个过程：

 
<BestCandidateSampler Declarations> =
    class BestCandidateSampler : public Sampler {
    public:
        < BestCandidateSampler Public Methods>
    private:
        < BestCandidateSampler Private Data>
    };

采样器把光栅空间中的当前像素位置存放在变量xTableCorner和yTableCorner中，而tableWidth是采样表所能覆盖的像素个数。变量tableOffset存放采样表中当前的偏置值，当该值前进到表的结尾，采样器将前进到图像的下一个区域。图7.33比较了最佳候选者模式和分层模式渲染棋盘的效果： (a)是分层模式，每个像素用了一个采样；(b)是最佳候选者模式，每个像素用了一个采样；(c)是分层模式，每个像素用了4个采样；(d)是最佳候选者模式，每个像素用了4个采样。虽然差别很细微，但仍可以看出在使用最佳候选者模式时，背景处的棋盘格的边走样化更轻微一些。还有，最佳候选者模式产生的噪声频率更高些，视觉上更容易接受。

 
图7.34显示了使用该模式对景深采样的效果。(a)是参考图像；(b)使用了低差异模式；(c)使用了最佳候选者模式。可以看出低差异模式的效果更佳，原因大约是在选择镜头采样时，是在一个固定区域中进行的，而根据实际所用的采样而定，该区域到图像上的映射可能远大于（或小于）单个像素面积。
 
<BestCandidateSampler Method Definitions> =
    BestCandidateSampler:: BestCandidateSampler(int xstart, int xend,
                             int ystart, int yend, int pixelSamples)
        : Sampler(xstart, xend, ystart, yend, pixelSamples) {
        tableWidth = (float) SQRT_SAMPLE_TABLE_SIZE / sqrtf(pixelSamples);
        xTableCorner = float(xPixelStart) - tableWidth;
        yTableCorner = float(yPixelStart);
        tableOffset = SAMPLE_TABLE_SIZE;
        < BestCandidateSampler constructor implementation>
    }

< BestCandidateSampler Private Data> =
    int tableOffset;
    float xTableCorner, yTableCorner, tableWidth;

下面把预先计算好的采样数据包括进来：

< BestCandidateSampler Private Data> +=
    static const float sampleTable[SAMPLE_TABLE_SIZE][5];
< BestCandidateSampler Method Definitions> +=
    #include "sampledata.cpp"

该采样器通常为积分器生成分层模式的采样（有一个例外，以后会解释）。实际上，这些模式对于正方形区域的分层效果更佳，所以RoundSize()将采样个数调整为一个整数的平方。

< BestCandidateSampler Public Methods> =
    int RoundSize(int size) const {
        int root = Ceil2Int(sqrtf((float)size - 0.5f));
        return root * root;
    }

BestCandidateSampler::GetNextSample()跟其它采样器中的方法基本相似，但不同的是，由于使用平铺方式，采样模式有时会伸出图像边界，这些出界的采样必须被忽略掉。

< BestCandidateSampler Method Definitions> +=
    bool BestCandidateSampler::GetNextSample(Sample *sample) {
    again:
        if (tableOffset = SAMPLE_TABLE_SIZE) {
            <Advance to next best-candidate sample table position>
        }
        <Compute raster sample from table>
        <Check sample against crop window, goto again if outside>
        <Compute integrator samples for best-candidate sample>
        ++tableOffset;
        return true;
    }

如果到达了采样表的尾部，采样器就试着向前移动xTableCorner (加上tableWidth)，如果出了图像边界，就要向前移动yTableCorner，如果在y方向上越过图像底端，就结束采样。

<Advance to next best-candidate sample table position> =
    tableOffset = 0;
    xTableCorner += tableWidth;
    if (xTableCorner >= xPixelEnd) {
        xTableCorner = float(xPixelStart);
        yTableCorner += tableWidth;
        if (yTableCorner >= yPixelEnd)
            return false;
    if (!oneSamples) {
        <Initialize sample tables and precompute strat2D values;
    }
    <Update sample shifts>
    <Generate SAMPLE_TABLE_SIZE-sized tables for signel samples>

这里在处理积分器的采样时，使用了一种混合方法：如果需要为每个图像采样生成单个的某种类型的采样，采样器就从一个被打乱次序的低差异序列中取一个值，所用的序列是一个为该段图像所生成的采样数组。如果需要生成多个采样，就要计算分层的采样。

使用这个方法有三个主要理由：首先，如果积分器为每个图像采样需要大量的采样，所使用最佳候选者采样表就会过大。第二，如果积分器采样的数目增加，使用分布并不那么好的采样对效果的影响也会减弱，因为积分器采样本身对当前的图像采样就已经有良好的覆盖了。第三，低差异序列跟某些像素重构滤波器一起使用时，会产生不寻常的人工痕迹。对于单个采样的情况而言，特别是PathIntegrator这样的需要许多单个采样的积分器，用低差异序列采样还是很值得的。

<BestCandidateSampler Private Data> +=
    float **oneDSamples, **twoDSamples;
    int *strat2D;

<BestCandidateSampler constructor implementation> =
    oneDSamples = twoDSamples = NULL;
    strat2D = NULL;

<Initialize sample tables and precompute strat2D values> =
    oneDSamples = new float * [sample->n1D.size()];
    for (u_int i = 0; i < sample->sample->n1D.size(); ++i) {
        oneDSamples = (sample->n1D == 1) ?
            new float[SAMPLE_TABLE_SIZE] : NULL;
    }
    twoDSamples = new float * [sample->n2D.size()];
    strat2D = new int [sample->n2D.size()];
    for (u_int i = 0; i < sample->sample->n2D.size(); ++i) {
        twoDSamples = (sample->n1D == 1) ?
            new float[2 * SAMPLE_TABLE_SIZE] : NULL;
        strat2D = Ceil2Int(sqrtf((float)sample->n2D - .5f));
    }

在单个采样情况下，低差异采样是用前面提到的LDShuffleScrambeld1D/2D工具函数生成的：

<Generate SAMPLE_TABLE_SIZE-sized tables for single samples> =
    for (u_int i = 0; i < sample->n1D.size(); ++i)
        if(sample->n1D == 1)
            LDShuffleScrambled1D(1, SAMPLE_TABLE_SIZE, oneDSamples);
    for (u_int i = 0; i < sample->n2D.size(); ++i)
        if(sample->n2D == 1)
            LDShuffleScrambled2D(1, SAMPLE_TABLE_SIZE, twoDSamples);

利用平铺式的采样模式有一个问题，就是有可能在模式的边缘产生细微的图像人工痕迹，这是因为相同的时间和镜头位置的采样在平铺时会重复出现。不仅如此，采样块的左上角总是有相同的时间和镜头值，等等。

一个解决方案是，在每次重用模式时对采样值进行变换。我们可以在每个维上使用Cranley-Patterson轮转变换：
        X'_i = ( Xi + ξ_i) mod 1

其中X_i是采样值，ξ_i是在0和1之间的随机值。因为我们所计算的采样模式具有环面拓扑结构，所以最终所得到的模式是具有良好分布和无缝的。时间和镜头的随机偏移量ξi在每次重用采样模式时都会被更新。

<Update sample shifts> =
    for (int i = 0; i < 3; i++)
        sampleOffset = RandomFloat();
<BestCandidateSampler Private Data> +=
    float sampleOffsets[3];

从采样表的位置来计算光栅空间采样位置需要一些简单的索引计算和比例计算。我们对图像采样并不使用Carnley-Patternson换位技术，因为这会在被重复的边界处产生不好的采样分布。保持图像采样的良好分布比减少关联更重要。我们在其它维上使用换位技术；WRAP宏定义保证结果落在0到1之间。

<Compute raster sample from table> =
    #define WRAP (x) ((x) > 1 ? ((x)-1) : (x))
    sample->imageX = xTableCorner + tableWidth * sampleTable[tableOffset][0];
    sample->imageY = yTableCorner + tableWidth * sampleTable[tableOffset][1];
    sample->time = WRAP(sampleOffset[0] + sampleTable[tableOffset][2]);
    sample->lensU = WRAP(sampleOffset[1] + sampleTable[tableOffset][3]);
    sample->lensV = WRAP(sampleOffset[2] + sampleTable[tableOffset][4]);

由于采样表可能超出图像平面的边界，所有有些采样会在采样区域之外，采样器需要处理这种情况：

<Check sample against crop window, goto again if outside> =
    if (         sample->imageX < xPixelStart ||
        sample->imageY >=  xPixelEnd ||
        sample->imageY< yPixelStart ||
        sample->imageY >=  yPixelEnd ) {
        ++tableOffset;
        goto again;
    }

前面解释过：对于积分器采样，如果需要一个采样，我们在一个预先计算好被扰乱的低差异序列中取值；否则就使用一个分层模式。

<Compute integrator samples for best-candidate sample> =
    for (u_int i = 0; i < sample->n1D.size(); ++i) {
        if(sample->n1D == 1)
            sample->oneD[0] = oneDSamples[tableOffset];
        else
            StratifiedSample1D(sampole->oneD, sample->n1D);
    for (u_int i = 0; i < sample->n2D.size(); ++i) {
        if(sample->n2D == 1) {
            sample->twoD[0] = oneDSamples[2 * tableOffset];
            sample->twoD[1] = oneDSamples[2 * tableOffset + 1];
        }
        else
            StratifiedSample2D(sampole->twoD, strat2D, strat2D);
    }

7.6 图像重构

有了精心挑选的图像采样，就要开发把采样及其辐射亮度值转换为供显示或存储的像素值的技术架构。根据信号处理理论，我们需要做下面的三件事：

1.    利用图像采样集合来重构图像采样L^~。
2.    对L^~做前置滤波，去掉超过Nyquist极限的频率。
3.    在像素位置上对L^~采样来计算出最后的像素值。

因为我们知道只在像素位置上对函数L~重新采样，所以就没有必要构造出函数的显式表达式出来，而是用一个滤波器将前两步合并。

回忆一下，如果用高于Nyquist频率的频率对原函数进行均匀采样并用sinc滤波器进行重构，所重构出的函数就可以跟原函数完全匹配 -- 这确实是件了不起的事情，因为我们所使用的仅仅是采样点而已。但由于图像函数几乎总是有更高的频率，以至于采样速率无法应付（如边缘上的情形，等等），我们不得不进行非均匀采样，把走样转化为噪声。

理想的重构是建立在均匀采样的基础之上的。虽然已经有些对非均匀采样理论进行扩展的尝试，但还没有令人接受的解决方法。还有，因为采样速率不足以覆盖整个函数，所以完美的重构也就不可能。最近的采样理论研究已经公认完美的重构在实际应用中通常是无法达到的，有了这个观点上的细微变化，就产生了关于重构的强有力的新技术（见Unser 2000 对这些技术的调查）。特别是研究重构理论的目标已经从完美重构转变到更现实的技术，也就是说，不管原函数是否有带宽局限，重构技术都要使重构函数跟原函数的差别最小化。

虽然pbrt所用的重构技术并不基于这些新技术，但是却也可以用来解释完美重构技术在图像合成过程中的不足之处。
为了重构像素值，我们来考虑对某个像素附近的采样进行插值的问题。为了计算像素的最终值I(x,y)，我们用下列的加权平均公式来进行计算：

 

其中L(xi,yi)是位于(xi,yi)的第i个采样的辐射亮度值，f是滤波函数。下图显示了位于(x,y)的像素，以及一个x,y方向上的范围分别是xWidth和yWidth的像素滤波器。所有位于滤波器范围之内的采样都会对像素值有所贡献，其贡献值有滤波函数值f(x-xi, y-yi)来决定。

 

这里sinc滤波器不是一个合适的选择：回忆一下，当函数有高过Nyquist极限的频率(Gibbs现象)时，理想的sinc滤波器容易产生振铃现象(Ringing)，也就是说图像中的边缘的附近会出现模糊的拷贝。还有，sinc滤波器有无限支撑(infinite support)：在有限距离之内，它不会衰减为0，因此对于每个像素，所有的图像采样都应该被滤波。实际上，并不存在最好的滤波函数。为特定的场景选择最佳的滤波器既需要量化的评估，又需要质量上的判断。

还有另一个问题可以影响对图像滤波器的选择：重构滤波器可能以某种令人诧异的方式对采样模式产生干扰。回忆一下LDSampler: 在单个像素的面积上，它可以生成分布非常良好的低差异模式，但像素的采样分布并没有考虑到临近像素的采样分布。当我们使用盒滤波器时，这个采样模式的效果非常出色，但当使用并非常数(constant)的滤波器并且跨越多个像素时，效果却不佳。这就产生左右为难的困境：我们最不想用的滤波器就是盒滤波器，但只有盒滤波器才能对LDSampler做出最佳的改进效果。基于对这些问题的考虑，pbrt提供了不同的滤波器插件。

7.6.1滤波函数

pbrt滤波器的实现都是从抽象类Filter来导出的，该抽象类提供了使用滤波函数f(x,y)的接口（见公式7.3）。类Film存放一个Filter指针，在把输出写入磁盘之前，用该滤波器对结果滤波。下图显示了不同滤波器的效果的比较：(a)是盒滤波器，(b)是Gaussian滤波器，(c)是Mitchell-Netravali滤波器。可以看出，Mitchell滤波后的图像最尖锐，Gaussian的图像最模糊，而盒滤波器的效果最差。

 

<Sampling Declarations > +=
    class Filter {
    public:
        <Filter Interface>
        <Filter Public Data>
    };

所有滤波器定义了一个宽度，超过该宽度值的函数值定义为0；该宽度在x和y方向上可以有不同的定义值。构造器把这些值及其倒数存放在相应的成员内，以备后用。滤波器在每个方向上的范围（即支撑）是相应宽度的两倍。见图：

 

<Filter Interface> =
    Filter(float xw, float yw)
        : xWidth(xw), yWidth(yw), invXWidth(1.f/xw), invYWidth(1.f/yw) {
    }

<Filter Public Data> =
    const float xWidth, yWidth;
    const float invXWidth, invYWidth;
Filter要实现的唯一函数是Evaluate()。它的变量x，y给出了采样点相对于滤波器中心的采样点位置。
<Filter Interface> +=
    virtual float Evaluate(float x, float y) const = 0;

盒滤波器

图形学中最常见的滤波器是盒滤波器，它对图像中的正方形区域中的所有采样都赋给相同的权值。虽然计算效率很高，但效果最糟糕。回忆一下第7.1节的讨论，盒滤波器可以让高频采样数据遗漏到重构值中。这就引起了后置走样(postaliasing) -- 虽然原采样值有足够高的频率来避免走样，但由于滤波效果很差而产生错误。
下图是盒滤波器的图形：

 

下图是利用盒滤波器对两个一维函数进行重构的结果。对于阶梯函数，盒滤波器做得很好，但对一个频率沿x轴增加正弦函数而言，则效果很差：不仅在频率低时重构的效果很差，而且当频率接近或超过Nyquist极限时效果更是糟糕透顶。

 

<Box Filter Declarations> =
    class BoxFilter : public Filter {
    public:
        BoxFilter(float xw, float yw) : Filter(xw, yw) { }
        float Evaluate(float x, float y) const;
    };

因为当(x,y)值超出滤波器范围时，我们不会调用Evaluate()函数，所以该函数返回1。

<Box Filter Method Definitions> =
    float BoxFilter::Evaluate(float x, float y) const {
        return 1.;
    }

三角形滤波器

三角形滤波器比盒滤波器的效果稍微好一些：在滤波器中心的采样有权值1，而其它的采样的权值沿中心两侧线性衰减。

 

<Triangle Filter Declarations>=
    class TriangleFilter : public Filter {
    public:
        TriangleFilter(float xw, float yw) : Filter(xw, yw) { }
        float Evaluate(float x, float y) const;
    };

<Triangle Filter Method Definitions> =
    float TriangleFilter::Evaluate(float x, float y) const {
        return max(0.f, xWidth - fabsf(x)) *
        max(0.f, yWidth - fabsf(y));
    }

高斯滤波器

与盒滤波器和三角形滤波器不同的是，高斯滤波器在实际应用中的效果相当不错。这个滤波器使用一个高斯凸起(Gaussian bump)，其中心位于像素的位置，并成放射式的对称分布。跟其它滤波器相比，高斯滤波器倾向于产生图像的轻微模糊，但是这种模糊实际上有助于去掉图像中剩余的走样。

 

<Gaussian Filter Declarations>=
    class GaussianFilter : public Filter {
    public:
        <GaussianFilter Public Methods>
    private:
        <GaussianFilter Private Data>
        <GaussianFilter Utility Functions>
    };

宽度为w的一维高斯滤波函数为：

 

其中α控制滤波器的衰减速率。α值越小，则衰减越慢，图像就越模糊。为了提供效率，构造器预先计算好了上面公式的在x,y方向上的常数项。

<GaussianFilter Public Methods> =
    GaussianFilter(float xw, float yw, float a)
    : Filter(xw, yw) {
    alpha = a;
    expX = expf(-alpha * xWidth * xWidth);
    expY = expf(-alpha * yWidth * yWidth);
}

<GaussianFilter Private Data>
    float alpha;
    float expX, expY;

因为2D高斯函数可以分离为两个1D高斯函数的积，所以下面的实现调用两次Gaussian()函数，并将结果相乘：

<Gaussian Filter Method Definitions> =
    float GaussianFilter::Evaluate(float x, float y) const {
    return Gaussian(x, expX) * Gaussian(y, expY);
}

<GaussianFilter Utility Functions> =
    float Gaussian(float d, float expv) const {
    return max(0.f, float(expf(-alpha * d * d) - expv));
}

Mitchell滤波器

设计滤波器是极其困难的，既需要数学上的分析，也需要感知上的经验。Mitchell和Netravali(1988)开发了一套参数化的滤波器函数来系统地探讨这个领域。通过分析试验者对不同参数值所产生的滤波图像的反应，他们开发出一个滤波器，它在两种常见的人工缺陷--即振铃现象和模糊现象之中做了很好的权衡。从下图的滤波器函数图形可以看出，有取负值的区域(称为negative lobes)。实际上这些负值区域改善了边缘的尖锐性，给出了更清楚的图像(减少了模糊性)。如果这些区域过大，则振铃现象就开始进入图像了。还有，虽然这可能产生负的图像值，但是这些负值最终要被转变为合法的范围中。

 

下图给出了滤波器对两个测试函数进行重构的效果。该滤波器对两者都表现得不错。对于(a)，振铃现象很微小，对正弦函数也做得很好，当然在函数频率增大的地方表现不佳，因为采样速率跟不上函数的变化。

 

<Mitchell Filter Declarations> =
    class MitchellFilter : public Filter {
    public:
        <MitchellFilter Public Methods>
    private:
        float B, C;
    };

Mitchell滤波器有两个参数: B和C。虽然两者可以取任意值，但Mitchell和Netravali建议它们保持 B + 2C = 1。

<MitchellFilter Public Methods> =
    MitchellFilter(float b, float c, float xw, float yw)
        : Filter(xw, yw) { B = b; C = c; }

跟其它许多2D图像滤波器一样（包括前面所介绍的高斯滤波器），Mitchell-Netravali滤波器是两个分别在x和y方向上的1D滤波器函数的乘积。

<Mitchell Filter Method Definitions> =
    float MitchellFilter::Evaluate(float x, float y) const {
        return Mitchell1D(x * invXWidth) * Mitchell1D(y * invYWidth);
    }

用于Mitchell滤波器的1D函数是定义在[-2, 2]上的偶函数。该函数用了在[0,1]上的三次多项式和[1,2] 上的三次多项式，并对平面x = 0进行镜像得到完整的函数。这两个多项式有B和C来控制，应保证在x=0, x= 1和x=2处的C0和C1的连续性。多项式定义如下：

 

<MitchellFilter Public Methods> =
    float Mitchell1D(float x) const {
        x = fabsf(2.f * x);
        if (x > 1.f)
            return ((-B - 6*C) * x*x*x + (6*B + 30*C) * x*x +
                (-12*B - 48*C) * x + (8*B + 24*C)) * (1.f/6.f);
        else
            return ((12 - 9*B - 6*C) * x*x*x +
                (-18 + 12*B + 6*C) * x*x +
                (6 - 2*B)) * (1.f/6.f);
    }

窗口Sinc滤波器

LanczosSincFilter类实现了基于sinc函数的滤波器。在实际应用中，sinc滤波器常常乘以一个随距离的增大而趋向于0的函数。这样滤波器函数就有了有限的范围，这对保证合理的性能还是很有必要的。另一个参数τ用来控制sinc函数的周期数。下图(a)表示了一个三个周期的sinc函数(实线)，和一个Lanczos开发的窗口函数(虚线)。该窗口函数定义如下：

 

两者的乘积就是下面(b)所表示的滤波器函数。

 

下面的图是对测试函数进行重构的结果。因为使用了窗口技术，跟用无限范围的sinc函数所重构出的函数相比，它所重构出的阶梯函数有更小的振铃现象。窗口sinc滤波器对正弦函数的重构也是极其成功的。

 

<Sinc Filter Declarations>=
    class LanczosSincFilter : public Filter {
    public:
    LanczosSincFilter(float xw, float yw, float t)
        : Filter(xw, yw) {
        tau = t;
    }

    float Evaluate(float x, float y) const;
    float Sinc1D(float x) const;
    private:
        float tau;
    };

跟其它滤波器一样，这个滤波器也是可分离的：

<Sinc Filter Method Definitions>=
    float LanczosSincFilter::Evaluate(float x, float y) const {
        return Sinc1D(x * invXWidth) * Sinc1D(y * invYWidth);
    }

下面的实现计算出sinc函数值，并跟Lanczos窗口函数相乘：

<Sinc Filter Method Definitions> +=
    float LanczosSincFilter::Sinc1D(float x) const {
        x = fabsf(x);
        if (x < 1e-5) return 1.f;
        if (x > 1.) return 0.f;
        x *= M_PI;
        float sinc = sinf(x * tau) / (x * tau);
        float lanczos = sinf(x) / x;
        return sinc * lanczos;
    }