Scape Terrain Editor - Giliam de

今天介绍的是Giliam de Carpentier的程序化地形方案，原文链接在文末有给出。

Scape是Carpentier为其学位论文而编写的地形编辑器，在这个编辑器中，允许用户使用程序化笔刷（Procedural Brushes，相对于传统的拉高压低等操作，程序化可以根据选择的地形类型的不同，生成不同细节的地形，比如山脉等）对地形进行编辑雕刻。

整个Terrain Editor系列包含了地形的渲染与地形数据的生成，渲染部分就是借助高度图加上一些传统的网格处理技术完成，没有太多值得介绍的内容，因此这里的主要内容就集中在地形数据的程序化生成算法上面，即如何通过程序化生成实现对自然地形的模拟。

1. 地形渲染

1.1 Shading

地形的着色是通过程序化的方式完成的，地形着色总共由四张贴图完成，通过高度以及地形梯度来控制不同位置的贴图混合权重，为了避免重复的纹样，这里还支持通过一张噪声贴图为不同的位置设置不同的混合方式，从而实现更为自然的混合效果。
为了解决根据顶点XZ坐标计算贴图采样UV坐标导致的拉伸，这里采用了经典的Tri-Planar采样方法

从而可以得到如下的效果：

具体的算法给出如下，大致思路就是使用法线平方作为混合权重，将三个方向轴采样之后的结果进行混合。

float3 uvwPos = uvwScaleFactor * worldPosition;
float3 weights = worldNormal * worldNormal;
float3 blendedColor = weights.xxx * tex2D(texSampler, uvwPos.yz).rgb +
                      weights.yyy * tex2D(texSampler, uvwPos.zx).rgb +
                      weights.zzz * tex2D(texSampler, uvwPos.xy).rgb;

解决投影扭曲的方案有多种，除了上面使用法线平方作为权重的方式，还有GPU Gem3的使用法线abs作为混合权重的方案，多种方案之间的效果对比可以参考如下的图片：

2. 地形生成

Scape的地形编辑器跟其他引擎的地形编辑一样，都是通过笔刷控制的，不同的是，这里的笔刷包含了内外两层，内层使用full strength，外层strength则从1.0过渡到0.0，此外这里还提供了不太常见的程序化笔刷功能，而程序化笔刷则是这篇文章的主要关注点。

Scape的程序化笔刷提供了众多的程序化算法，但是这些算法都是基于2D Perlin噪声的，这是因为Perlin噪声具有如下的一些优点：

• 噪声是band-limited的
• 计算方便快捷
• 可以实现局部以及并行计算（即无过多关联数据需要考虑）

2.1 Perlin噪声

Perlin噪声支持任意维度的输入与输出，对于地形生成，我们只需要传入一个2D坐标p，外加一个1D的噪声种子，就能够收获一个对应位置的1D输出，其算法给出如下：

uniform sampler2D samplerPerlinPerm2D;
uniform sampler2D samplerPerlinGrad2D;

float perlinNoise(float2 p, float seed)
{
    // Calculate 2D integer coordinates i and fraction p.
    float2 i = floor(p);
    float2 f = p - i;

    // Get weights from the coordinate fraction
    float2 w = f * f * f * (f * (f * 6 - 15) + 10);
    float4 w4 = float4(1, w.x, w.y, w.x * w.y);

    // Get the four randomly permutated indices from the noise lattice nearest to
    // p and offset these numbers with the seed number.
    float4 perm = tex2D(samplerPerlinPerm2D, i / 256) + seed;

    // Permutate the four offseted indices again and get the 2D gradient for each
    // of the four permutated coordinates-seed pairs.
    float4 g1 = tex2D(samplerPerlinGrad2D, perm.xy) * 2 - 1;
    float4 g2 = tex2D(samplerPerlinGrad2D, perm.zw) * 2 - 1;

    // Evaluate the four lattice gradients at p
    float a = dot(g1.xy, f);
    float b = dot(g2.xy, f + float2(-1,  0));
    float c = dot(g1.zw, f + float2( 0, -1));
    float d = dot(g2.zw, f + float2(-1, -1));

    // Bi-linearly blend between the gradients, using w4 as blend factors.
    float4 grads = float4(a, b - a, c - a, a - b - c + d);
    float n = dot(grads, w4);

    // Return the noise value, roughly normalized in the range [-1, 1]
    return n * 1.5;
}

samplerPerlinPerm2D用于对一个精心设计的4通道彩色贴图（分辨率为256x256，如下图所示）进行采样。这张贴图是一张LUT，对于一个给定的整数坐标i，会返回四个坐标位置(i.x, i.y), of (i.x, i.y + 1), of (i.x + 1, i.y) and of (i.x + 1, i.y + 1) 的hash值，这四个hash值就分别存储于四个通道中。 这张贴图采样模式是点采样，而非双线性采样，在采样的过程中需要关闭mipmap。

Permutation Texture

samplerPerlinGrad2D 的使用逻辑类似，不同的是对于某个采样点，会将上一个贴图的采样结果根据seed对其整数坐标偏移，从而得到两个相邻的整数坐标，用这两个坐标作为输入分别对gradient贴图进行一次lookup，得到两个梯度向量。对于包含四个顶点的cell来说，我们只需要两个梯度向量就能够计算出对应cell的梯度. 根据这些数据就可以计算p点的高度，最后各点的高度会使用quintic插值算法实现插值并归一化。跟perm贴图一样，gradient贴图也是需要repeating+point sampling + no mipmapping。

Gradient Texture

2.2 FBM

FBM (fractional Brownian motion)是对Perlin噪声进行组合的最简单算法，一些更为复杂的turbulence算法都是以这个算法为基础完成的。其基本原理就是对不同频率的噪声以不同的振幅进行累加，这种算法在很多方面都有应用，比如在体积云实现中，用于对Voronoi噪声进行处理以实现更为自然的效果。

FBM累加中的每一个不同频率的噪声被称为一个Octave，不同的Octave负责不同尺寸的feature，算法实现给出如下：

float turbulence(float2 p, float seed, int octaves, float lacunarity = 2.0, float gain = 0.5)
{
  float sum = 0;
  float freq = 1.0, amp = 1.0;
  for (int i=0; i < octaves; i++)
  {
    float n = perlinNoise(p*freq, seed + i / 256.0);
    sum += n*amp;
    freq *= lacunarity;
    amp *= gain;
  }
  return sum;
}

算法中的seed通常是一个[0, 1]范围内的浮点数，而为了避免重复的pattern，lacunarity数值通常会被设定为一个低于2.0的浮点数（比如1.92），为了能够覆盖基本上所有的feature，octave数目通常会少于log(terrain_width) / log(lacunarity)，对于一个1024x1024尺寸的地形而言，我们通常就需要大约10个octave的perlin噪声累加；gain数值是不同octave噪声之间的振幅差异，用于控制不同的粗糙粒度。

2.3 Billowy turbulence

如果用下面的billowedNoise替换之前的PerlinNoise的函数实现，我们就可以得到一些凸起部分具有billow（鼓鼓的）效果而凹陷部分带有腐蚀Erosion效果且具有较为尖锐褶痕（crease）的地形结果，这个噪声在很多文献中就直接被称为turbulence function，实际上就是Perlin噪声的绝对值。

float billowedNoise(float2 p, float seed)
{
    return abs(perlinNoise(p, seed));
}

2.4 Ridged turbulence

通过对Turbulence Function取反，就可以得到山脉一样尖锐的山脊效果

float ridgedNoise(float2 p, float seed)
{
    return 1.0f-abs(perlinNoise(p, seed));
}

2.5 IQ Noise

上面的算法是通过对Perlin噪声进行累加得到的，而这种做法的弊端在于波峰跟波谷的尺寸都是近似的，而这跟自然中的地形风貌不太一致。

解决这个问题的一种思路是，在对多个octave的噪声进行累加前，先根据粗糙的octave噪声的输出来确定更为精细的octave噪声的振幅，从而得到更为异构的地形效果，这个算法最开始是由I. Quillez提出的：

float iqTurbulence(float2 p, float seed, int octaves, float lacunarity = 2.0, float gain = 0.5)
{
    float sum = 0.5;
    float freq = 1.0, amp = 1.0;
    float2 dsum = float2(0,0);
    for (int i=0; i < octaves; i++)
    {
        float3 n = perlinNoisePseudoDeriv(p*freq, seed + i / 256.0);
        dsum += n.yz;
        // the higher the octave, the smaller the amplitude
        sum += amp * n.x / (1 + dot(dsum, dsum));
        freq *= lacunarity;
        amp *= gain;
    }
    return sum;
}

上述算法中，n是一个通过perlinNoisePseudoDeriv 函数计算得到的三维向量，其中n.x就是此前算法中输出的Perlin噪声，也就是此前Turbulence Function中的n。而n.yz则有点类似于当前采样点下Perlin噪声的局部微分，这里说有点类似，是因为原作者在推导微分计算公式的时候存在错误，但是输出的效果却歪打正着，正好是我们所需要的效果。

在octave迭代的过程中，这个微分数值会累加起来对振幅进行抑制，从而减少斜坡（高微分值）区域的细节凹凸效果。下面给出perlinNoisePseudoDeriv 算法的具体实现：

float3 perlinNoisePseudoDeriv(float2 p, float seed)
{
    // Calculate 2D integer coordinates i and fraction p.
    float2 i = floor(p);
    float2 f = p - i;

    // Get weights from the coordinate fraction
    float2 w = f * f * f * (f * (f * 6 - 15) + 10);
    float4 w4 = float4(1, w.x, w.y, w.x * w.y);

    // Get pseudo derivative weights
    float2 dw = f * f * (f * (30 * f - 60) + 30);

    // Get the four randomly permutated indices from the noise lattice nearest to
    // p and offset these numbers with the seed number.
    float4 perm = tex2D(samplerPerlinPerm2D, i / 256) + seed;

    // Permutate the four offseted indices again and get the 2D gradient for each
    // of the four permutated coordinates-seed pairs.
    float4 g1 = tex2D(samplerPerlinGrad2D, perm.xy) * 2 - 1;
    float4 g2 = tex2D(samplerPerlinGrad2D, perm.zw) * 2 - 1;

    // Evaluate the four lattice gradients at p
    float a = dot(g1.xy, f);
    float b = dot(g2.xy, f + float2(-1,  0));
    float c = dot(g1.zw, f + float2( 0, -1));
    float d = dot(g2.zw, f + float2(-1, -1));

    // Bi-linearly blend between the gradients, using w4 as blend factors.
    float4 grads = float4(a, b - a, c - a, a - b - c + d);
    float n = dot(grads, w4);

    // Calculate pseudo derivates
    float dx = dw.x * (grads.y + grads.w*w.y);
    float dy = dw.y * (grads.z + grads.w*w.x);

    // Return the noise value, roughly normalized in the range [-1, 1]
    // Also return the pseudo dn/dx and dn/dy, scaled by the same factor
    return float3(n, dx, dy) * 1.5;
}

2.6 Swiss turbulence

下面两张图中，Procedural erosion给出的是通过gidged turbulence函数生成的地形效果，这种地形看起来并不太差，但是缺少了地形成型的那种历史演变信息。通常来说，自然地形中会在地形的斜坡上生成槽沟，而在山谷中的地形会较为平缓，而这些效果在这个地形中均没有，虽然我们可以通过对heightfield进行处理来模拟这种侵蚀效果，但是需要经过大量迭代，效率较低，幸运的是，我们可以通过算法进行fake，下面Procedural ridged turbulence就给出了对应的输出效果。

Procedural erosion Procedural ridged turbulence

用来实现这个效果的算法叫做SwissTurbulence，这个算法依然是通过将多个Perlin噪声组合起来实现的，具体算法逻辑给出如下：

float swissTurbulence(float2 p, float seed, int octaves,
                      float lacunarity = 2.0, float gain = 0.5,
                      float warp = 0.15>)
{
     float sum = 0;
     float freq = 1.0, amp = 1.0;
     float2 dsum = float2(0,0);
     for(int i=0; i < octaves; i++)
     {
         float3 n = perlinNoiseDeriv((p + warp * dsum)*freq, seed + i);
         sum += amp * (1 - abs(n.x));
         // specific modifications
         dsum += amp * n.yz * -n.x;
         freq *= lacunarity;
         amp *= gain * saturate(sum);
    }
    return sum;
}

可以看到，相对于此前的Ridged Turbulence函数，这里只是增加了一个dsum的计算逻辑，这部分后面会介绍。perlinNoiseDeriv 函数跟perlinNoisePseudoDeriv 函数很相似，不同的是，这里的yz分量返回是的正确推导的微分结果，而x分量还是Perlin噪声。

整个算法结构还是跟之前一样，通过对噪声应用FBM来完成，不同的是，由于这里希望生成带有Ridge效果的地形，因此这里的累加是针对1-abs(n.x)而非n.x的，通过这种方式可以生成基础的ridged turbulence效果。

Step 1: Basic ridged turbulence

dsum跟sum一样，都是一个加权求和项，不同的是sum是针对噪声的求和，而dsum则是针对微分数值的求和。最终这个2D的求和项会被更精一级octave用作噪声offset，这个offset将被更精一级octave中的采样点用来查找最近的ridge数据，从而导致如下图所示的在ridge在斜坡上的拉长效果。这个拉长效果的幅度受warp参数的影响，需要注意的是，按照chain rule（复合求导）来说，1-abs(n.x)的梯度公式应该写成n.xy-sign(n.x) 而非n.yz-n.x，不过为了避免sign函数的不连续导致的异常，这里就直接使用后者进行模拟了。

Step 2: Distort by gradient

前面说过，山谷地形应该要比较平缓，为了实现这个效果，这里的amp计算方式做了轻微修正，amp不再是只乘上每个octave的gain，还会乘上一个临时的sum变量（即上面算法中的saturate(sum)），而这个做法会使得更精细的octave的细节数据在地形高度降低的时候fade out，而在地形高度较高区域则继续维持甚至更为突出。而为了对整体振幅缩小的补偿，这里使用的gain参数应该要比正常情况下要稍大，比如从0.5改成0.6。

Step 3: Smoothen valleys

另外，为了使ridge以及slope效果具有更为随机自然的效果，这里还可以对输入的p做一次distortion。这里的做法是将p分别传入到两个不同的turbulence function（跟这里的turbulence function不一样的两个函数），将输出的结果组合成一个新的p。这里对p进行distortion扰动的函数所使用的gain，octave，scale等参数都不必跟SwissTurbulence函数中一致，且为了得到较好的效果，这里扰动函数中的octave数目不能太高（通常4就够了）。

Step 4: Distort by noise

下面给出perlinNoiseDeriv函数的实施逻辑。

float3 perlinNoiseDeriv(float2 p, float seed)
{
    // Calculate 2D integer coordinates i and fraction p.
    float2 i = floor(p);
    float2 f = p - i;

    // Get weights from the coordinate fraction
    float2 w = f * f * f * (f * (f * 6 - 15) + 10); // 6f^5 - 15f^4 + 10f^3
    float4 w4 = float4(1, w.x, w.y, w.x * w.y);

    // Get the derivative dw/df
    float2 dw = f * f * (f * (f * 30 - 60) + 30); // 30f^4 - 60f^3 + 30f^2

    // Get the derivative d(w*f)/df
    float2 dwp = f * f * f * (f * (f * 36 - 75) + 40); // 36f^5 - 75f^4 + 40f^3

    // Get the four randomly permutated indices from the noise lattice nearest to
    // p and offset these numbers with the seed number.
    float4 perm = tex2D(samplerPerlinPerm2D, i / 256) + seed;

    // Permutate the four offseted indices again and get the 2D gradient for each
    // of the four permutated coordinates-seed pairs.
    float4 g1 = tex2D(samplerPerlinGrad2D, perm.xy) * 2 - 1;
    float4 g2 = tex2D(samplerPerlinGrad2D, perm.zw) * 2 - 1;

    // Evaluate the four lattice gradients at p
    float a = dot(g1.xy, f);
    float b = dot(g2.xy, f + float2(-1,  0));
    float c = dot(g1.zw, f + float2( 0, -1));
    float d = dot(g2.zw, f + float2(-1, -1));

    // Bi-linearly blend between the gradients, using w4 as blend factors.
    float4 grads = float4(a, b - a, c - a, a - b - c + d);
    float n = dot(grads, w4);

    // Calculate the derivatives dn/dx and dn/dy
    float dx = (g1.x + (g1.z-g1.x)*w.y) + ((g2.y-g1.y)*f.y - g2.x +
               ((g1.y-g2.y-g1.w+g2.w)*f.y + g2.x + g1.w - g2.z - g2.w)*w.y)*
               dw.x + ((g2.x-g1.x) + (g1.x-g2.x-g1.z+g2.z)*w.y)*dwp.x;
    float dy = (g1.y + (g2.y-g1.y)*w.x) + ((g1.z-g1.x)*f.x - g1.w + ((g1.x-
               g2.x-g1.z+g2.z)*f.x + g2.x + g1.w - g2.z - g2.w)*w.x)*dw.y +
               ((g1.w-g1.y) + (g1.y-g2.y-g1.w+g2.w)*w.x)*dwp.y;

    // Return the noise value, roughly normalized in the range [-1, 1]
    // Also return the pseudo dn/dx and dn/dy, scaled by the same factor
    return float3(n, dx, dy) * 1.5;
}

2.7 Jordan turbulence

下面两张图给出的高度图效果依然是通过对Perlin噪声进行混合输出的，没有使用任何积分方式的迭代算法。

跟SwissTurbulence函数一样，这里也同样使用了perlinNoiseDeriv 函数对octave进行扰动跟扭曲，但是却输出了一个完全不同的地形效果：

float jordanTurbulence(float2 p, float seed, int octaves, float lacunarity = 2.0,
                       float gain1 = 0.8, float gain = 0.5,
                       float warp0 = 0.4, float warp = 0.35,
                       float damp0 = 1.0, float damp = 0.8,
                       float damp_scale = 1.0)
{
    float3 n = perlinNoiseDeriv(p, seed);
    float3 n2 = n * n.x;
    float sum = n2.x;
    float2 dsum_warp = warp0*n2.yz;
    float2 dsum_damp = damp0*n2.yz;

    float amp = gain1;
    float freq = lacunarity;
    float damped_amp = amp * gain;

    for(int i=1; i < octaves; i++)
    {
        n = perlinNoiseDeriv(p * freq + dsum_warp.xy, seed + i / 256.0);
        n2 = n * n.x;
        sum += damped_amp * n2.x;
        dsum_warp += warp * n2.yz;
        dsum_damp += damp * n2.yz;
        freq *= lacunarity;
        amp *= gain;
        damped_amp = amp * (1-damp_scale/(1+dot(dsum_damp,dsum_damp)));
    }
    return sum;
}

跟其他扰动算法不同的是，这里最粗糙一级的octave噪声是在for循环之外计算的。这种做法的好处是最粗糙一级octave的振幅将跟其他octave的噪声独立开来，当然更好的一种控制方式是将每一级octave的权重都分开处理，通过一个数组进行控制，从而可以节省GPU cycle，同时提供更为灵活的控制。

首先，最终输出的是噪声之和，这里的噪声之和是n2.x，也就是(n.x)^2，通过这个平方，使得输出的噪声具有billowy的效果，而这也是JordanTurbulence跟SwissTurbulence不同的原因之一了。当warp/warp0/damp/damp0都是0的时候，这时候dsum_warp/dsum_damp也都是0，这时候噪声效果就相当于普通Turbulence噪声的平方，效果就跟下面一样。

Step 1: Squared noise turbulence

根据复合求导算法，(n.x)^2的梯度是2.0n.xn.yz或2.0*n2.yz，因此dsum_damp会与叠加上一个scale因子的梯度进行累加，而因为上一级的噪声结果比较平缓，因此这个结果会用于对下一级的octave的振幅进行平滑处理，这个效果看起来就像是与thermal erosion的组合。

跟SwissTurbulence不同，JordanTurbulence根据梯度而非高度来对下一级噪声的振幅进行dampen（平滑处理），这也就意味着不论是山峰还是山谷，平缓区域处的噪声细节都会比较少。其中damp0控制着最粗糙一级噪声对后续噪声damping效果的强度，而damp则控制着当前一级噪声对后面噪声的damping效果强度。被damping的振幅大小是受damp_scale控制的。

Step 2: Gradient-based damping

dsum_warp变量的计算方式跟dsum_damp类似，不过其使用目的不一样。这个变量的作用是用于创建类似下斜槽沟的下过，用于模拟fluvial erosion，因此其控制参数需要保持跟dsum_damp不同。作为副作用，这个变量还会导致一个对所有feature的一个squash（压挤）效果。

warp0用于控制最粗糙一级的噪声对所有其他级数噪声的作用，而warp则控制着上一级噪声对余下所有噪声的作用。

Step 3: Gradient-based warping

这里的所有步骤的操作方式跟SwissTurbulence差不多，不过最终输出了完全不同的地形效果。前面说过SwissTurbulence中会需要对输入的p进行扰动处理，这里也是一样的，同样需要通过两个不同Turbulence函数对p进行扰动并组合成新的p来计算后续效果，其中一个可选的Turbulence函数可以使用前面的SwissTurbulence，这个函数可以缓解最终JordanTurbulence中波峰距离比较接近的瑕疵。

Step 4: Adding 2D noise to p

2.8 笔刷编辑管线

Scape的编辑管线主要包含三个阶段：

1. 捕捉用户输入，并将之转化为笔刷操作存入到一个操作队列中

• 用户的输入经过投影后被翻译成cubic spline上的一系列采样点
• 根据用户配置，对cubic spline以一定的间隔进行采样
• 上一步中的采样点会被当成笔刷操作实例添加到队列中

2. 对笔刷操作实例进行处理

• 每个笔刷操作实例被分割成一个或者多个paged笔刷实例
• 通过一个scheduler对待处理的paged笔刷实例batch进行处理
• 对page的处理会将对应的地形tile设置为invalidation
• 当某次编辑对应的所有笔刷操作实例都处理完成后，老的地形就会被压缩存储，同时更新一下undo stack

3. 更新terrain geometry

• 使用新的heightfield page作为输入，完成之前设置为invalidation的tile的重新生成

Page是一个具有固定尺寸的方块，Page跟Tile不是一个概念，最大的区别是Page是没有LOD概念的，通过Page可以避免地形整体尺寸过大导致的操作的不方便（如Undo/Redo等的不方便）。

参考资料

[1]. Scape - Rendering the terrain

[2]. Scape - Procedural basics

[3]. Scape - Procedural extensions

[4]. Scape - Brush pipeline

[5]. Scape - Overview and downloads

[6]. Effective GPU-based synthesis and editing of realistic heightfields - Giliam de Carpentier