重要性采样(Importance Sampling)简介和简单样例实现

在渲染领域，重要性采样这个术语是很常见的，但它究竟是什么呢？我们首先考虑这样的一种情况：

如果场景里有一点P，我们想计算P点的最终颜色，根据全局照明的概念，P点的颜色是由所有投射到P点的所有光线所影响的。但是我们在利用蒙特卡洛积分去计算投射光线的总和会出现一个问题，如上图场景里有两盏灯，当蒙特卡洛采样较少的情况，可能导致灯光方向未被采样，这样计算的结果会和实际产生巨大的偏差，尤其是在采样函数(f(X))有突峰的情况。如下图所示：

而在采样函数是常数时，则完全没有该问题：

当然，我们几乎不会使用是常数的采样函数。可是我们可以这样计算蒙特卡洛积分：

我们寻找一个pdf(X)函数是f(X)的倍数，这样采样函数会变成如下图所示的常数了：

该图正好是我们上面两盏灯场景的光线强度采样函数，当我们取 ${f}'(X)$ 为 ${f}(X)$ 的倍数，则我们的采样就能像采样常数那样进行。当我们的pdf选择上图中所示的函数时，我们的X则正好有在两盏灯光附近采样较多样本，其他地方采样较少的样本的特点。这就是我们所说的重要性采样。虽然这不是均匀采样，但是我们在 ${f}(X)$ 大的时候，我们的pdf值也大， ${f}(X)$ 小的时候，pdf值也小，类似相互补足了，也导致我们的函数最终能够正确的收敛，而且具有更小的方差。但实际上很难找到完全是倍数的pdf，但我们可以找形状十分类似的pdf。

我们现在用一个例子：

我们使用蒙特卡洛积分和重要性采样来计算该积分。虽然该积分可以很快速的用计算获得：

但是我们的目的是为了验证重要性采样的快速收敛性。我们选择两个pdf进行对比，如下图所示：

一个是 $p(x)=1/\pi$ 的均匀分布，另一个则是 ${p}'(x)=8x/\pi^{2}$ 形状与 $f(x)$ 类似的分布。在我们用代码验证之前，我们需要考虑的是我们通常用的drand48()是返回一个均匀分布的随机值，如何获得想要的分布的随机值呢？

我们先从正态分布的例子：

假设我们有一组数据，火者到站的时间与预估的到站时间差符合正态分布。我们这么利用drand48()获得正态分布的随机值呢？我们可以想到正态分布的CDF函数，因为其单调递增且值域为[0,1]：

我们的均匀随机分布正好可以是其CDF的y值，我们用y取计算X，就可以获得我们需要的正态分布(可以是任意)的随机值了。利用y计算X有两种方法，一种是直接利用公式反推的方法，另一种首先用程序根据pdf计算CDF表(离散形式)，然后利用drand48函数获取随机值y并在CDF表中索引并计算出X。示例代码如下：

#include <cstdlib> 
#include <cstdio> 
#include <cstring> 
#include <cmath> 
#include <algorithm> 
 
// 正态分布pdf
float pdf(const float &x) 
{ return 1 / sqrtf(2 * M_PI) * exp(-x * x * 0.5); } 
 
float sample(float *cdf, const uint32_t &nbins, const float &minBound, const float &maxBound) 
{ 
    float r = drand48(); 
    float *ptr = std::lower_bound(cdf, cdf + nbins + 1, r); 
    int off = std::max(0, (int)(ptr - cdf - 1)); 
    float t = (r - cdf[off]) / (cdf[off + 1] - cdf[off]); 
    float x = (off + t) / (float)(nbins); 
    return minBound + (maxBound - minBound) * x; 
} 
 
int main(int argc, char ** argv) 
{ 
     srand48(13); 
     // 生成CDF表
     int nbins = 32; 
     float minBound = -5, maxBound = 5; 
     float cdf[nbins + 1], dx = (maxBound - minBound) / nbins, sum = 0; 
     cdf[0] = 0.f; 
     for (int n = 1; n < nbins; ++n) { 
         float x = minBound + (maxBound - minBound) * (n / (float)(nbins)); 
         float pdf_x = pdf(x) * dx; 
         cdf[n] = cdf[n - 1] + pdf_x; 
         sum += pdf_x; 
     } 
     cdf[nbins] = 1; 
 
    // 我们的模拟
    int numSims = 100000; 
    int numBins = 100; // to collect data on our sim 
    int bins[numBins]; // to collect data on our sim 
    memset(bins, 0x0, sizeof(int) * numBins); // 设定初始值0 
    const float dist = 10; // 10 km 
    for (int i = 0; i < numSims; ++i) { 
	   float diff = sample(cdf, nbins, minBound, maxBound); // 利用y获取X
	   int whichBin = (int)(numBins * (diff - minBound) / (maxBound - minBound)); 
	   bins[whichBin]++; 
	   float time = 30 + diff; 
        float speed = 60 * dist / time; 
    } 
    sum = 0; 
    for (int i = 0; i < numBins; ++i) { 
        float r = bins[i] / (float)numSims; 
        printf("%f %f\n", 5 * (2 * (i /(float)(numBins)) -1), r); 
        sum += r; 
    } 
    fprintf(stderr, "sum %f\n", sum); 
    return 0; 
}

本文不采用上述代码中的方法(上述代码的方法只是作为一个实例)，而是直接使用公式反推。我们的公式：

所以我们可以用drand48()获取 $\xi$ 的值，然后计算x，这样x就是满足pdf为 ${p}'(x)=8x/\pi^{2}$ 分布的随机值了。我们把使用了重要性采样和使用均匀分布的积分作对比：

#include <cstdio> 
#include <cstdlib> 
#include <cmath> 
int main(int argc, char **argv) 
{ 
    srand48(13); 
    int N = 16; 
    for (int n = 0; n < 10; ++n) { 
        float sumUniform = 0, sumImportance = 0; 
        for (int i = 0; i < N; ++i) { 
            float r = drand48(); 
            sumUniform += sin(r * M_PI * 0.5); 
            float xi = sqrtf(r) * M_PI * 0.5; // this is our X_i 
            sumImportance += sin(xi) / ((8 * xi) / (M_PI * M_PI)); 
        } 
        sumUniform *= (M_PI * 0.5) / N; 
        sumImportance *= 1.f / N; 
        printf("%f %f\n", sumUniform, sumImportance); 
    } 
    return 0; 
}

我们输出的结果如下表：