实时体积云渲染(地平线)：二.Perlin噪声和Worley噪声

Perlin噪声

Perlin噪声同样是网格点噪声的一种，不同于之间在网格点生成随机值的白噪声，Perlin噪声在网格点生成一个随机的单位化三维向量，我们可以直接生成3个随机[0,1]的浮点数，然后重新映射到[-1,1]即可。代码如下：

unsigned seed = 2;
std::mt19937 generator(seed);
std::uniform_real_distribution distribution;
auto dice = std::bind(distribution, generator);
float gradientLen2;
for (unsigned i = 0; i < tableSize; ++i) { gradients[i] = Vec3f(2 * dice() - 1, 2 * dice() - 1, 2 * dice() - 1); gradientLen2 = gradients[i].length2(); gradients[i] /= sqrtf(gradientLen2); // normalize gradient permutationTable[i] = i;
}

由于噪声函数返回的是实数，这里的计算方式就是利用网格点上的梯度和网格点到需要计算噪声的位置P的向量做点乘，以3D噪声为例，我们将利用P点和周围8个网格点梯度向量做点乘，然后利用三线性插值计算噪声值。代码如下：

class PerlinNoise
{ static const unsigned tableSize = 256; static const unsigned tableSizeMask = tableSize - 1; Vec3f gradients[tableSize]; unsigned permutationTable[tableSize * 2]; PerlinNoise() { unsigned seed = 2016; std::mt19937 generator(seed); std::uniform_real_distribution distribution; auto dice = std::bind(distribution, generator); float gradientLen2; for (unsigned i = 0; i < tableSize; ++i) { do { gradients[i] = Vec3f(2 * dice() - 1, 2 * dice() - 1, 2 * dice() - 1); gradientLen2 = gradients[i].length2(); } while (gradientLen2 > 1); gradients[i] /= sqrtf(gradientLen2); // normalize gradient permutationTable[i] = i; } std::uniform_int_distribution distributionInt; auto diceInt = std::bind(distributionInt, generator); // 创建 permutation 表for (unsigned i = 0; i < tableSize; ++i) std::swap(permutationTable[i], permutationTable[diceInt() & tableSizeMask]); // 扩大表的范围至 [256:512]for (unsigned i = 0; i < tableSize; ++i) { permutationTable[tableSize + i] = permutationTable[i]; } } virtual ~PerlinNoise() {} /* inline */ int hash(const int &x, const int &y, const int &z) const { return permutationTable[permutationTable[permutationTable[x] + y] + z]; } public: float eval(const Vec3f &p) const { int xi0 = ((int)std::floor(p.x)) & tableSizeMask; int yi0 = ((int)std::floor(p.y)) & tableSizeMask; int zi0 = ((int)std::floor(p.z)) & tableSizeMask; int xi1 = (xi0 + 1) & tableSizeMask; int yi1 = (yi0 + 1) & tableSizeMask; int zi1 = (zi0 + 1) & tableSizeMask; float tx = p.x - ((int)std::floor(p.x)); float ty = p.y - ((int)std::floor(p.y)); float tz = p.z - ((int)std::floor(p.z)); float u = smoothstep(tx); float v = smoothstep(ty); float w = smoothstep(tz); // 获取周围8个网格点的梯度const Vec3f &c000 = gradients[hash(xi0, yi0, zi0)]; const Vec3f &c100 = gradients[hash(xi1, yi0, zi0)]; const Vec3f &c010 = gradients[hash(xi0, yi1, zi0)]; const Vec3f &c110 = gradients[hash(xi1, yi1, zi0)]; const Vec3f &c001 = gradients[hash(xi0, yi0, zi1)]; const Vec3f &c101 = gradients[hash(xi1, yi0, zi1)]; const Vec3f &c011 = gradients[hash(xi0, yi1, zi1)]; const Vec3f &c111 = gradients[hash(xi1, yi1, zi1)]; // 生成指向p的向量float x0 = tx, x1 = tx - 1; float y0 = ty, y1 = ty - 1; float z0 = tz, z1 = tz - 1; Vec3f p000 = Vec3f(x0, y0, z0); Vec3f p100 = Vec3f(x1, y0, z0); Vec3f p010 = Vec3f(x0, y1, z0); Vec3f p110 = Vec3f(x1, y1, z0); Vec3f p001 = Vec3f(x0, y0, z1); Vec3f p101 = Vec3f(x1, y0, z1); Vec3f p011 = Vec3f(x0, y1, z1); Vec3f p111 = Vec3f(x1, y1, z1); // 线性插值float a = lerp(dot(c000, p000), dot(c100, p100), u); float b = lerp(dot(c010, p010), dot(c110, p110), u); float c = lerp(dot(c001, p001), dot(c101, p101), u); float d = lerp(dot(c011, p011), dot(c111, p111), u); float e = lerp(a, b, v); float f = lerp(c, d, v); return lerp(e, f, w);}

值得注意的是，这里我们使用随机表大小为256。这代表我们的噪声以256为一个周期，即[-256,0]，[256,512]都与[0,256]一致。2D的周期显示结果如下图：

我们这里还使用了Permutation方法。

Permutation

如果要生成一维的噪声，我们设定表大小为256，那么2维的表该设定多大呢？256*256？其实这里只需要256+256的大小，而3维也只需要256+256。这就是用到了Permutation方法。

为了解决这一限制，我们要做的就是确保噪声函数的任何输入值（浮点，2D或3D点）都可以映射到该随机值数组中的一个且只有一个值。如您所知，每个输入点与我们的噪声函数的坐标都转换为整数。该整数用于查找随机值数组。肯·佩林（Ken Perlin）将数字0到255存储在另一个由256个整数组成的数组中，并对这个数组的条目进行混洗（置换），随机地交换它们中的每一个。实际上，Perlin将此数组的大小扩展到512，并将索引范围0到255的值复制到索引256到511。换句话说，将数组的前半部分复制到后半部分。具体代码为：

        // 生成 permutation 表for (unsigned i = 0; i < tableSize; ++i) std::swap(permutationTable[i], permutationTable[diceInt() & tableSizeMask]); // 扩大表的范围至 [256:512]for (unsigned i = 0; i < tableSize; ++i) { permutationTable[tableSize + i] = permutationTable[i]; }

由于我们的x，y，z都对255求“与”运算，所以我们的hash函数为：

    int hash(const int &x, const int &y, const int &z) const { return permutationTable[permutationTable[permutationTable[x] + y] + z]; }

这样便可以始终用[0,512]作为索引，获取新的[0,512]的随机值。

噪声值计算

计算噪声值的方式非常简单，以2D为例子，就是用最近的4个网格点到点P的向量和他们自身的梯度求点积，所得再使用双线性插值计算(代码中为3D，所以用三线性插值)。为了使网格点影响在近处更大，我们对双线性插值计算之前，进行smoothstep重映射。

我们运行后的结果如下图：

FBM

在使用噪声时，我们经常会看到FBM这个词眼，它的全名为分型布朗运动(Fractal Brownian Motion)，其特点就是函数处处连续且不可导。在噪声中，我们可以理解为无数不同频的噪声叠加，我们常用方式是波高每次迭代*0.5，频率每次迭代*2。对上图迭代10次，有如下效果：

代码如下：

                for(int fbmi=0;fbmi<FBMnum;fbmi++){//计算出的噪声值为[-1,1]重映射为[0,1]noiseMap[indexM]+=(FBMh*(pnose.eval(Vec3f(i,k,j)*(FBMf / 64.0))+1.0)*0.5);FBMf*=2;FBMh/=2;}//累加过程可能会超过1，因此需要重新映射为[0,1]if(noiseMap[indexM]>fbmMax){fbmMax = noiseMap[indexM];}for (uint32_t k = 0; k <height*width; ++k) {noiseMap[k] /= fbmMax;}

均匀分布的梯度

我们上述生成梯度的方法其实并不能获得均匀分布的梯度。我们在正方体的每个面上均匀撒点会导致正方体角上团聚大量的样本。如果我们使用ϕ是[0,2π]，θ是[0，π]，并写出代码：

float phi = 2 * drand48() * M_PI;
float theta = drand48() * M_PI;
float x = cos(phi) * sin(theta);
float y = sin(phi) * sin(theta);
float z = cos(theta);

然而这会产生如下左图情况：

我们需要了解如何才能真正的生成球面均匀分布的样本。

为了获得球面上均匀采样的样本，我们有：

其中p(w)为概率密度函数，因为是均匀的，所以每个样本的概率一致，所以p(w)=C，C为常数。解得：

我们用ϕ是[0,2π]，θ是[0，π]来表示p，有：

根据立体角公式：

带入：

可得：

可以利用求解边缘密度的方法求解p(ϕ)和p(θ)。

由于ϕ和θ相互独立，所以p(ϕ)可以直接求值：

现在我们有了ϕ和θ的PDF，即概率密度函数，我们可以求解它们的CDF，即概率积累函数了，p(θ)的CDF为：

求解变限积分：

同理，p(ϕ)的CDF为：

我们就获得了p(ϕ)和p(θ)的CDF，我们对其求逆，便可以使用均匀的[0,1]的分布获取ϕ和θ分布了。θ求逆过程如下：

ϕ求逆过程如下：

至此，我们就能获得球面均匀分布的采样了，更改噪声生成代码：

for (unsigned i = 0; i < tableSize; ++i) { float theta = acos(2 * dice() - 1); float phi = 2 * dice() * M_PI; float x = cos(phi) * sin(theta); float y = sin(phi) * sin(theta); float z = cos(theta); gradients[i] = Vec3f(x, y, z); permutationTable[i] = i;
}

Worley噪声

worley噪声在网格点设置的噪声值不同，它是给区域划分网格，在每个网格内部随机生成一个种子点，我们任意点的噪声值通过和周围最近的种子点计算得出。示例图如下图：

小白点为生成的随机种子点，我们发现距离种子点越近，则噪声值越小。根据该规则，我们代码为：

class WorleyNoise
{
public:WorleyNoise(const unsigned &seed = 2016){std::mt19937 generator(seed);std::uniform_real_distribution<float> distribution;auto dice = std::bind(distribution, generator);for (unsigned i = 0; i < tableSize; ++i) {float x = dice();float y = dice();float z = dice();positions[i] = Vec3f(x, y, z);permutationTable[i] = i;}std::uniform_int_distribution<unsigned> distributionInt;auto diceInt = std::bind(distributionInt, generator);// create permutation tablefor (unsigned i = 0; i < tableSize; ++i)std::swap(permutationTable[i], permutationTable[diceInt() & tableSizeMask]);// extend the permutation table in the index range [256:512]for (unsigned i = 0; i < tableSize; ++i) {permutationTable[tableSize + i] = permutationTable[i];}}virtual ~WorleyNoise() {}float eval(const Vec3f &p) const{int xi0 = ((int)std::floor(p.x)) & tableSizeMask;int yi0 = ((int)std::floor(p.y)) & tableSizeMask;int zi0 = ((int)std::floor(p.z)) & tableSizeMask;float m_dist=1.0;for(int i = -1; i < 2; i++){for (int j = -1; j < 2; j++) {for(int k = -1; k < 2; k++){Vec3i neighborI = Vec3i((xi0+i)& tableSizeMask,(yi0+j)& tableSizeMask,(zi0+k)& tableSizeMask);Vec3i neighborPos=Vec3i(((int)std::floor(p.x))+i,((int)std::floor(p.y))+j,((int)std::floor(p.z))+k);Vec3f neighborP=Vec3f(neighborPos.x,neighborPos.y,neighborPos.z)+positions[hash(neighborI.x, neighborI.y, neighborI.z)];float dist = distance(p,neighborP);m_dist=fmin(dist, m_dist);}}}return m_dist;}private:/* inline */uint8_t hash(const int &x, const int &y, const int &z) const{return permutationTable[permutationTable[permutationTable[x] + y] + z];}static const unsigned tableSize = 128;static const unsigned tableSizeMask = tableSize - 1;Vec3f positions[tableSize];unsigned permutationTable[tableSize * 2];
};

由于这里我们设置的噪声表大小为128，因此我们的周期即为128，我们为了保证在边界周围计算距离时保证原来距离，所以我们使用neighborP去保存原始位置。完成效果如下图：

我们在生成云使用的往往是1-worley噪声，即返回值为1-m_dist。我们修改后的效果为：

使用FBM后，生成的效果：

Perlin-Worley噪声

Perlin-Worley噪声就是两者的结合，我们使用公式：

clamp(remap(PerlinNoise, worleyNoise, 1.0, 0.0, 1.0));

就能得到Perlin-Worley噪声了。效果如下：

生成和使用3D噪声

我们如何生成3D的噪声并在着色器中使用呢？我们可以将3D的噪声图平铺在2D的纹理中，形式如下：

纹理格式为(s,t,r)。随着r增加，平铺方式为从左到右，从上倒下。以生成Perlin-Worley噪声为例，我们代码为：

    WorleyNoise noise(7);PerlinNoise pnose(2);// output noise map to PPMconst uint32_t width = 128, height = 128, depth = 8;float *noiseMap = new float[width * height * depth * depth];float *worleyMap= new float[width * height * depth * depth];int FBMnum=8;float fbmMax=0,worMax=0;for (uint32_t k = 0; k < depth*depth; ++k) {for (uint32_t j = 0; j < height; ++j) {for (uint32_t i = 0; i < width; ++i) {float FBMf=1.0,FBMh=1.0;//FBM频率，波高uint32_t indexM=(k/depth)*width*height*depth+width*depth*j;indexM+=(k%depth)*width+i;for(int fbmi=0;fbmi<FBMnum;fbmi++){worleyMap[indexM]+=FBMh*noise.eval(Vec3f(i,k,j)*(FBMf / 64.0));noiseMap[indexM]+=(FBMh*(pnose.eval(Vec3f(i,k,j)*(FBMf / 64.0))+1.0)*0.5);FBMf*=2;FBMh/=2;}if(noiseMap[indexM]>fbmMax){fbmMax = noiseMap[indexM];}if(worleyMap[indexM]>worMax){worMax= worleyMap[indexM];}}}}for (uint32_t k = 0; k < depth*depth*height*width; ++k) {noiseMap[k] /= fbmMax;worleyMap[k]/=worMax;noiseMap[k]=clamp(remap(noiseMap[k], worleyMap[k], 1.0, 0.0, 1.0));}std::ofstream ofs;ofs.open("./noise2.ppm", std::ios::out | std::ios::binary);ofs << "P6\n" << width*depth << " " << height*depth << "\n255\n";for (unsigned k = 0; k < depth*depth*height*width; ++k) {unsigned char n = static_cast<unsigned char>(noiseMap[k] * 255);ofs << n << n << n;}ofs.close();delete[] noiseMap;

生成的噪声效果如下：

我们接下来需要用glTexImage3D绑定纹理，绑定时，OpenGL规定3D顺序为x->y->z，因此纹理同样为s->t->r。所以我们纹理像素的读取顺序为下图：

生成texture的代码如下：

GLuint general3DTex(const std::string& name,unsigned int sizeT){GLuint textures;int width,height,nrChannels;glGenTextures(1, &textures);glBindTexture(GL_TEXTURE_3D, textures);glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_REPEAT);glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_REPEAT);glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_REPEAT);glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);unsigned char* image=stbi_load(name.c_str(), &width, &height, &nrChannels, 0);unsigned char* fimage=new unsigned char[width*height*nrChannels];//云噪声图为128int tex2Dsize=sizeT;int depth=(width/tex2Dsize)*(width/tex2Dsize);int tex2Dblock=width/tex2Dsize;for (int i=0; i<width; i++) {for (int j=0; j<height; j++) {unsigned int indexM=j/tex2Dsize*tex2Dblock+i/tex2Dsize;indexM=indexM*tex2Dsize*tex2Dsize+j%tex2Dsize*tex2Dsize+i%tex2Dsize;for(int k=0;k<nrChannels;k++){fimage[nrChannels*indexM+k]=image[nrChannels*(width*j+i)+k];}}}if(!image){std::cout<<"texture error to load"<<std::endl;return -1;}if(nrChannels==3){glTexImage3D(GL_TEXTURE_3D, 0, GL_RGB, tex2Dsize, tex2Dsize,depth,0, GL_RGB, GL_UNSIGNED_BYTE, fimage);}else if (nrChannels==4){glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA, tex2Dsize, tex2Dsize,depth, 0, GL_RGBA, GL_UNSIGNED_BYTE, fimage);}stbi_image_free(image);delete []fimage;glBindTexture(GL_TEXTURE_2D, 0);return textures;
}

我们可视化生成的3DPerlin-Worley噪声如下图所示：