ShaderToy常见代码解释 (收录中)

The Book Of Shaders 的 smoothstep 函数绘制 Shader 解析

#ifdef GL_ES
precision mediump float;
#endif

#define PI 3.14159265359

float plot(vec2 st, float pct) {
  return smoothstep(pct-0.02, pct, st.y) -
         smoothstep(pct, pct+0.02, st.y);
}

void main() {
  vec2 st = gl_FragCoord.xy / iResolution;

  // smooth interpolation between 0.1 and 0.9
  float y = smoothstep(0.1, 0.9, st.x);

  vec3 color = vec3(y);

  float pct = plot(st, y);
  color = (1.0 - pct) * color + pct * vec3(0.0, 1.0, 0.0);

  gl_FragColor = vec4(color, 1.0);
}

float plot(vec2 st, float pct) {
  return abs(pct - st.y) < 0.01 ? 1.0: 0.0;
  /* 或者使用 step 函数
     step(abs(pct - st.y), 0.01);

     type = float/vec2/vec3/vec4

     type step(type edge, type x)

     equals

     returnValue[i] = edge[i] > x[i] ? 0.0: 1.0;
   */
  // return step(abs(pct - st.y), 0.01) ? 1.0: 0.0;
}

ShaderToy 的编程思路总结

高维随机函数

float rand(vec2 st) {
  return fract(sin(dot(st.xy, vec2(12.9898, 78.233))) * 43758.5453123);
}

向量值函数形式的随机函数

// https://www.shadertoy.com/view/lldyDn
vec2 rand(vec2 p)
{
  vec3 p3 = fract(vec3(p.xyx) * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yzx+19.19);
  return fract((p3.xx+p3.yz) * p3.zy);
}

vec3 rand(vec3 v) {
  vec3 p = vec3(dot(v, vec3(127.1, 311.7, 74.7)), dot(v, vec3(269.5, 183.3, 246.1)), dot(v, vec3(113.5, 271.9, 124.6)));
  return fract(sin(p) * 43758.5453123);
}

vec3 rand(vec2 v) {
  v -= 0.00001727738;
  float z = fract(tan(v.x*3.14+v.y))*10.751135;
  vec3 a = fract(vec3(v.x/1.000383,v.y/1.000401,z)*3.14)*10.8372611876;
  float w = fract(dot(a, vec3(1.3774, 8.2446, 2.2883)))*10.0973326;
  return fract(7.817082*tan(9.726749*reflect(a, vec3(w, z, fract(z*43.28)*10.114865))));
}

三角函数的缺陷

// https://www.shadertoy.com/view/llGSzw
float hash11 ( uint n ) {
  // integer hash copied from Hugo Elias
  n = (n << 13U) ^ n;
  n = n * (n * n * 15731U + 789221U) + 1376312589U;
  return float( n & uint(0x7fffffffU))/float(0x7fffffff);
}

vec3 hash13( uint n )
{
  // integer hash copied from Hugo Elias
  n = (n << 13U) ^ n;
  n = n * (n * n * 15731U + 789221U) + 1376312589U;
  uvec3 k = n * uvec3(n,n*16807U,n*48271U);
  return vec3( k & uvec3(0x7fffffffU))/float(0x7fffffff);
}

// https://www.shadertoy.com/view/4tXyWN
uint hash21( uvec2 p )
{
  p *= uvec2(73333,7777);
  p ^= (uvec2(3333777777)>>(p>>28));
  uint n = p.x*p.y;
  return float(n^(n>>15)) * (1.0 / float(0xffffffffU));
}

// https://www.shadertoy.com/view/XlXcW4
vec3 hash33( uvec3 p ) {
  const uint k = 1103515245U;    // GLIB C
  //const uint k = 134775813U;   // Delphi and Turbo Pascal
  //const uint k = 20170906U;    // Today's date (use three days ago's dateif you want a prime)
  //const uint k = 1664525U;     // Numerical Recipes
  p = ((p>>8U)^p.yzx) * k;
  p = ((p>>8U)^p.yzx) * k;
  p = ((x>>8U)^p.yzx) * k;

  return vec3(x) * (1.0 / float(0xffffffffU));
}

//  1 out, 1 in...
float chash11(float p)
{
  p = fract(p * .1031);
  p *= p + 33.33;
  p *= p + p;
  return fract(p);
}

//----------------------------------------------------------------------------------------
//  1 out, 2 in...
float chash12(vec2 p)
{
  vec3 p3  = fract(vec3(p.xyx) * .1031);
  p3 += dot(p3, p3.yzx + 33.33);
  return fract((p3.x + p3.y) * p3.z);
}

//----------------------------------------------------------------------------------------
//  1 out, 3 in...
float chash13(vec3 p3)
{
  p3  = fract(p3 * .1031);
  p3 += dot(p3, p3.zyx + 31.32);
  return fract((p3.x + p3.y) * p3.z);
}

//----------------------------------------------------------------------------------------
//  2 out, 1 in...
vec2 chash21(float p)
{
  vec3 p3 = fract( p * vec3(.1031, .1030, .0973) );
  p3 += dot(p3, p3.yzx + 33.33);
  return fract((p3.xx+p3.yz)*p3.zy);

}

//----------------------------------------------------------------------------------------
///  2 out, 2 in...
vec2 chash22(vec2 p)
{
  vec3 p3 = fract(vec3(p.xyx) * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yzx+33.33);
  return fract((p3.xx+p3.yz)*p3.zy);

}

//----------------------------------------------------------------------------------------
///  2 out, 3 in...
vec2 chash23(vec3 p3)
{
  p3 = fract(p3 * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yzx+33.33);
  return fract((p3.xx+p3.yz)*p3.zy);
}

//----------------------------------------------------------------------------------------
//  3 out, 1 in...
vec3 chash31(float p)
{
  vec3 p3 = fract(p * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yzx+33.33);
  return fract((p3.xxy+p3.yzz)*p3.zyx);
}

//----------------------------------------------------------------------------------------
///  3 out, 2 in...
vec3 chash32(vec2 p)
{
  vec3 p3 = fract(vec3(p.xyx) * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yxz+33.33);
  return fract((p3.xxy+p3.yzz)*p3.zyx);
}

//----------------------------------------------------------------------------------------
///  3 out, 3 in...
vec3 chash33(vec3 p3)
{
  p3 = fract(p3 * vec3(.1031, .1030, .0973));
  p3 += dot(p3, p3.yxz+33.33);
  return fract((p3.xxy + p3.yxx)*p3.zyx);

}

//----------------------------------------------------------------------------------------
// 4 out, 1 in...
vec4 chash41(float p)
{
  vec4 p4 = fract( p * vec4(.1031, .1030, .0973, .1099));
  p4 += dot(p4, p4.wzxy+33.33);
  return fract((p4.xxyz+p4.yzzw)*p4.zywx);

}

//----------------------------------------------------------------------------------------
// 4 out, 2 in...
vec4 chash42(vec2 p)
{
  vec4 p4 = fract(vec4(p.xyxy) * vec4(.1031, .1030, .0973, .1099));
  p4 += dot(p4, p4.wzxy+33.33);
  return fract((p4.xxyz+p4.yzzw)*p4.zywx);

}

//----------------------------------------------------------------------------------------
// 4 out, 3 in...
vec4 chash43(vec3 p)
{
  vec4 p4 = fract(vec4(p.xyzx)  * vec4(.1031, .1030, .0973, .1099));
  p4 += dot(p4, p4.wzxy+33.33);
  return fract((p4.xxyz+p4.yzzw)*p4.zywx);
}

//----------------------------------------------------------------------------------------
// 4 out, 4 in...
vec4 chash44(vec4 p4)
{
  p4 = fract(p4  * vec4(.1031, .1030, .0973, .1099));
  p4 += dot(p4, p4.wzxy+33.33);
  return fract((p4.xxyz+p4.yzzw)*p4.zywx);
}

float mod289(float x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec2 mod289(vec2 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec3 mod289(vec3 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

vec4 mod289(vec4 x) {
  return x - floor(x * (1.0 / 289.0)) * 289.0;
}

// mod(34.0 * x * x + 10.0 * x, 289.0)
/* 34.0 * x * x + 10.0 * x 为非线性函数, 该函数具备雪崩效应(avalanche effect),
   即输入发生微小改变会导致输出发生巨大且不可预测的变化, 或者说具备较大的非常数导数 */
float permute(float x) {
     return mod289(((x*34.0)+10.0)*x);
}

vec2 permute(vec2 x) {
     return mod289(((x*34.0)+10.0)*x);
}

vec3 permute(vec3 x) {
     return mod289(((x*34.0)+10.0)*x);
}

vec4 permute(vec4 x) {
     return mod289(((x*34.0)+10.0)*x);
}

值噪声

// 来自 https://thebookofshaders.com/11/
float rand(vec2 st) {
  return fract(sin(dot(st.xy, vec2(12.9898, 78.233))) * 43758.5453123);
}

float noise (in vec2 st) {
  vec2 i = floor(st);
  vec2 f = fract(st);

  // Four corners in 2D of a tile which is created according to st.
  float a = rand(i);
  float b = rand(i + vec2(1.0, 0.0));
  float c = rand(i + vec2(0.0, 1.0));
  float d = rand(i + vec2(1.0, 1.0));

  // Smooth Interpolation

  // Cubic Hermine Curve.  Same as SmoothStep()
  vec2 u = f*f*(3.0-2.0*f);
  // u = smoothstep(0.,1.,f);

  // Mix 4 coorners percentages
  return mix(a, b, u.x) +
          (c - a)* u.y * (1.0 - u.x) +
          (d - b) * u.x * u.y;
}

void main() {
  vec2 st = gl_FragCoord.xy/iResolution.xy;

  // Scale the coordinate system to see
  // some noise in action
  vec2 pos = vec2(st * 5.0);

  // Use the noise function
  float n = noise(pos);

  gl_FragColor = vec4(vec3(n), 1.0);
}

梯度噪声

// 2D 向量值函数形式的随机函数
vec2 rand(vec2 st){
    st = vec2(dot(st, vec2(127.1,311.7)),
              dot(st, vec2(269.5,183.3)));
    return -1.0 + 2.0 * fract(sin(st) * 43758.5453123);
    // 归一化
}

// 类似前面 2D 噪声的生成, 但是 4 个角的随机值需要一番点积运算
float noise(vec2 st) {
  vec2 i = floor(st);
  vec2 f = fract(st);
  // vec2 u = f * f * (3.0 - 2.0 * f);
  vec2 u = smoothstep(vec2(0.0), vec2(1.0), f);

  // a,b,c,d 是四个网格上的梯度
  // va,vb,vc,vd 是四个角指向
  vec2 a = rand(i);
  vec2 va = f - vec2(0.0, 0.0);
  vec2 b = rand(i + vec2(1.0, 0.0));
  vec2 vb = f - vec2(1.0, 0.0);
  vec2 c = rand(i + vec2(0.0, 1.0));
  vec2 vc = f - vec2(0.0, 1.0);
  vec2 d = rand(i + vec2(1.0, 1.0));
  vec2 vd = f - vec2(1.0, 1.0);

  float dotA = dot(a, va);
  float dotB = dot(b, vb);
  float dotC = dot(c, vc);
  float dotD = dot(d, vd);

  float ab = mix(dotA, dotB, u.x);
  float cd = mix(dotC, dotD, u.x);

  return mix(ab, cd, u.y);
}

float rand(float p) {
  return fract(sin(p) * 43758.5453123);
}

float noise(float p) {
  float i = floor(p);
  float f = fract(p);
  float u = smoothstep(0.0, 1.0, f);

  // 这里的 rand(i) * 2.0 - 1.0 是要对梯度 rand(i) 进行归一化使其范围在 [-1, 1] 之间
  float a = rand(i) * 2.0 - 1.0;
  float va = f - 0.0;
  float b = rand(i + 1.0) * 2.0 - 1.0;
  float vb = f - 1.0;

  float dotA = a * va;
  float dotB = b * vb;

  return 2.4 * mix(dotA, dotB, u);
}

单形噪声

• 单形
  

  在可视化空间时, 会用选择一种几何图形来重复地对空间进行填充, 或说, 选择一种几何图形对空间进行细分(tesselation), 从而形成坐标网格.

  为了方便说明, 这里把正交空间上的网格叫做正交网格, 单形空间上的网格叫做单形网格.

  在细分一些空间时, 几何图形的选择会有很多种, 而几何结构最简单的那一种选择就叫做单形(simplex shape), 也可以理解为空间的单位.

  一维空间的单形是等长的线段(一维空间只有这种选择).

  对于二维空间, 通常会用正方形对它进行细分, 但正方形本身就能够划分成两个三角形.

  因此, 三角形才是二维空间的单形, 二维空间的最优单形是等边三角形, 这意味着要把正方形沿着主对角线 \(y = x\) 错切成菱形, 使得划分所得的三角形等边, 2 个等边三角形可以组合成一个菱形.

  

  Figure 20: 等边三角形作为单形, 填充二维空间

  对于三维空间, 它的单形是四面体(tetrahedron), 正式来说应该是正四面体, 可通过把对立方体沿着主对角线 \(x = y = z\) 进行错切, 再根据对角线划分得到 6 个正四面体.

  

  Figure 21: 正四面体作为单形, 填充三维空间

  对于四维空间, 它的单形很难可视化, 不过它有 5 个角, 并且 24 个单形可以组合成一个错切的四维超立方体(hypercube).

  总的来说, \(N\) 维空间的单形有 \(N + 1\) 个角, \(N!\) 个单形可以填充完一个有着 \(2^{N}\) 个角的错切 \(N\) 维超立方体.

  在传统噪声算法中, 随着走向更高维度, 为每个角计算随机梯度是一个复杂度为 \(O(2^{N})\) 的问题,

  如果是对单形的每个角度计算随机梯度, 那么计算复杂度会变成 \(O(N^2)\), 这就是单形的好处.

• 用求和替代插值
  

  随着维度走高, 传统噪声函数除了计算复杂度变高外, 噪声函数的解析导数也会越来越难求解, 噪声函数的解析导数有很多用处:

  包括凹凸贴图(bump mapping), 位移贴图(displacement mapping), 解析抗锯齿(analytical antialiasing)以及使用旋度噪声(curl nosie)的粒子动画.

  为了解决这个问题, 单形噪声不再像梯度噪声那样根据两角的贡献值进行插值, 而是对每个角的贡献值进行求和作为噪声值.

  

  在单形噪声中, 角的贡献值的计算过程如下:

  \([\max(0, r^2 - |\vec{d}|^2)]^4 \times \vec{d} \cdot \vec{g}\), 其中 \(\vec{d}\) 是角指向, \(\vec{g}\) 是梯度值, \(\max(0, r^2 - |\vec{d}|^2)\) 是衰减函数, \(|\vec{d}|^2\) 是角指向的模长.

  衰减函数让贡献值随着点与角之间的距离增大而减少, 通常 \(r^2 = 0.5\), 也就是以角为圆心且半径为 \(\sqrt{0.5} \approx 0.7\) 的圆,

  一般来说不会让这个圆覆盖同一个单形里的其它单形角, 因此, 半径 \(r\) 应该设定为等边三角形的高, 而不是等边三角形的边长.

  在后面的实现二维单形噪声时会详细说明为什么是 \(r^2 = 0.5\).

  

• 决定点所处的单形, 并且决定角的遍历顺序
  

  在传统噪声算法中需要判断点处于哪个网格中, 在单形噪声算法中则需要判断点处于哪个单形中.

  二维空间的单行是由菱形划分的等边三角形, 分上三角形和下三角形, 判断点在那个单形也就是判断它在上三角形还是下三角形.

  

  在二维空间上, 如果点 \((x ,y)\) 满足 \(x \gt y\), 那么点处在下三角形中; 反之, 点处在上三角形中.

  如果点处于下三角形中, 那么它的单形角的遍历顺序为 \((0, 0)\), \((1, 0)\) 和 \((1, 1)\).

  如果点处于上三角形中, 那么它的单形角的遍历顺序为 \((0, 0)\), \((0, 1)\) 和 \((1, 1)\).

  决定这个顺序很简单, 为了说明, 这里采用 \((i, j)\) 来作为单形的相对角坐标:

  要对点 \((x, y)\) 的分量进行从大到小排序, 这个顺序决定了角相对坐标分量加 1 的顺序.

  如果点的分量顺序是 \(x \gt y\), 那么是对上一个角相对坐标 \(a\) 的 \(i\) 分量先加 1, 得到角相对坐标 \(b\), 然后相对 \(b\) 的 \(j\) 分量加 1, 得到角相对坐标 \(c\):

  \(a = (i, j), b = (i + 1, j), c = (i + 1, j + 1)\).

  如果点的分量顺序是 \(y \gt z\), 那么是对上一个角相对坐标 \(a\) 的 \(j\) 分量先加 1, 得到角相对坐标 \(b\), 然后相对 \(b\) 的 \(i\) 分量加 1, 得到角相对坐标 \(c\):

  \(a = (i, j), b = (i, j + 1), c = (i + 1, j + 1)\).

  对于更高空间的点 \((x, y, \dots)\) 和角相对坐标 \((x, y, \dots)\) 同理.

  下图是三维空间中判断所处单形的方法.

  

• 2D 单形噪声实现
  

  这里是解读 Simplex Noise Demystified 原文中 2D Simplex Noise 的实现:

  // 2D simplex noise
  public static double noise(double xin, double yin) {
      double n0, n1, n2; // Noise contributions from the three corners
      // Skew the input space to determine which simplex cell we're in
      final double F2 = 0.5*(Math.sqrt(3.0)-1.0);
      double s = (xin+yin)*F2; // Hairy factor for 2D
      int i = fastfloor(xin+s);
      int j = fastfloor(yin+s);
  
      final double G2 = (3.0-Math.sqrt(3.0))/6.0;
      double t = (i+j)*G2;
      double X0 = i-t; // Unskew the cell origin back to (x,y) space
      double Y0 = j-t;
      double x0 = xin-X0; // The x,y distances from the cell origin
      double y0 = yin-Y0;
  
      // For the 2D case, the simplex shape is an equilateral triangle.
      // Determine which simplex we are in.
      int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
      if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1)
      else {i1=0; j1=1;}
      // upper triangle, YX order: (0,0)->(0,1)->(1,1)
      // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
      // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
      // c = (3-sqrt(3))/6
  
      double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords
      double y1 = y0 - j1 + G2;
      double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords
      double y2 = y0 - 1.0 + 2.0 * G2;
  
      // Work out the hashed gradient indices of the three simplex corners
      int ii = i & 255;
      int jj = j & 255;
      int gi0 = perm[ii+perm[jj]] % 12;
      int gi1 = perm[ii+i1+perm[jj+j1]] % 12;
      int gi2 = perm[ii+1+perm[jj+1]] % 12;
  
      // Calculate the contribution from the three corners
      double t0 = 0.5 - x0*x0-y0*y0;
      if(t0<0) n0 = 0.0;
      else {
          t0 *= t0;
          n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient
      }
      double t1 = 0.5 - x1*x1-y1*y1;
      if(t1<0) n1 = 0.0;
      else {
          t1 *= t1;
          n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
      }double t2 = 0.5 - x2*x2-y2*y2;
      if(t2<0) n2 = 0.0;
      else {
          t2 *= t2;
          n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
      }
      // Add contributions from each corner to get the final noise value.
      // The result is scaled to return values in the interval [-1,1].
      return 70.0 * (n0 + n1 + n2);
  }
  

  这段代码远没有作者口中的那么可读, 很多细节都没有说明, 不过确实如他所说的很方便的移植到其它语言上.

  不过也不能怪作者, 因为原本 Ken Perlin 的实现就很难理解了.

  • 把正交网格变换成单形网格
    

    在实现单形噪声时, 可以 自由 选择使用任何单形网格, 这里选择了 Simplex Noise Demystified 里的单形网格, 如下图:

    

    这个单形网格是这么得到的:

    把 \(x\) 轴 \(\left( \begin{array}{c} 1 \\ 0 \end{array} \right)\) 往逆时针方向旋转 \(15^{\circ}\) 得到 \(u\) 轴, 再把 \(y\) 轴 \(\left( \begin{array}{c} 0 \\ 1 \end{array} \right)\) 往顺时针方向旋转 \(15^{\circ}\) 得到 \(v\) 轴, 使得 \(u\) 和 \(v\) 两轴之间的夹角为 \(60^{\circ}\); 最后对 \(u\) 和 \(v\) 两轴进行统一伸缩.

    这个结论是通过 Simplex Noise Demystified 以下代码反推的, 原文没有明确说明如何得到的单形网格, 因此, 结论可能有误.

    final double F2 = 0.5*(Math.sqrt(3.0)-1.0);
    double s = (xin+yin)*F2; // Hairy factor for 2D
    int i = fastfloor(xin+s);
    int j = fastfloor(yin+s);
    
    final double G2 = (3.0-Math.sqrt(3.0))/6.0;
    double t = (i+j)*G2;
    double X0 = i-t; // Unskew the cell origin back to (x,y) space
    double Y0 = j-t;
    double x0 = xin-X0; // The x,y distances from the cell origin
    double y0 = yin-Y0;
    

    为了方便讨论, 假设 \(\vec{u} = \left( \begin{array}{c} a \\ c \end{array} \right)\) 和 \(\vec{v} = \left( \begin{array}{c} b \\ d \end{array} \right)\) 就是旋转得到的 \(u\) 和 \(v\) 轴, 并且都是单位向量.

    接下来计算得出 \(u\) 和 \(v\) 的矩阵:

    \(\cos \frac{\pi}{12} = \cos(\frac{\pi}{4} - \frac{\pi}{6}) = \cos\frac{\pi}{4} \cos\frac{\pi}{6} + \sin\frac{\pi}{4} \sin\frac{\pi}{6} = \frac{\sqrt{3}}{2 \sqrt{2}} + \frac{1}{2 \sqrt{2}} = \frac{1}{2 \sqrt{2}}(\sqrt{3} + 1)\)

    \(\cos (-\frac{\pi}{12}) = \cos(\frac{\pi}{6} - \frac{\pi}{4}) = \cos\frac{\pi}{6}\cos\frac{\pi}{4} + \sin\frac{\pi}{6}\sin\frac{\pi}{4} = \frac{\sqrt{3}}{2 \sqrt{2}} + \frac{1}{2 \sqrt{2}} = \frac{1}{2 \sqrt{2}}(\sqrt{3} + 1)\)

    \(\sin \frac{\pi}{12} = \sin(\frac{\pi}{4} - \frac{\pi}{6}) = \sin\frac{\pi}{4} \cos\frac{\pi}{6} - \cos\frac{\pi}{4} \sin\frac{\pi}{6} = \frac{\sqrt{3}}{2 \sqrt{2}} - \frac{1}{2 \sqrt{2}} = \frac{1}{2 \sqrt{2}}(\sqrt{3} - 1)\)

    \(\sin (-\frac{\pi}{12}) = \sin(\frac{\pi}{6} - \frac{\pi}{4}) = \sin\frac{\pi}{6} \cos\frac{\pi}{4} - \cos\frac{\pi}{6} \sin\frac{\pi}{4} = \frac{1}{2 \sqrt{2}} - \frac{\sqrt{3}}{2 \sqrt{2}} = \frac{1}{2 \sqrt{2}}(1 - \sqrt{3})\)

    \(R(\frac{\pi}{12}) = \left( \begin{array}{c} \cos(\frac{\pi}{12}) & -\sin(\frac{\pi}{12}) \\ \sin(\frac{\pi}{12}) & \cos(\frac{\pi}{12}) \end{array} \right) = \frac{1}{2 \sqrt{2}} \left( \begin{array}{c} \sqrt{3} + 1 & 1 - \sqrt{3} \\ \sqrt{3} - 1 & \sqrt{3} + 1 \end{array} \right)\)

    \(R(-\frac{\pi}{12}) = R(\frac{\pi}{12})^{-1} = R(\frac{\pi}{12})^{T}\)

    \(\vec{x}\) 轴变换得到 \(\vec{u}\) 轴: \(\vec{u} = R(\frac{\pi}{12}) \left( \begin{array}{c} 1 \\ 0 \end{array} \right) = \frac{1}{2\sqrt{2}}\left( \begin{array}{c} \sqrt{3} + 1 \\ \sqrt{3} - 1 \end{array} \right)\).

    \(\vec{y}\) 轴变换得到 \(\vec{v}\) 轴: \(\vec{v} = R(-\frac{\pi}{12}) \left( \begin{array}{c} 0 \\ 1 \end{array} \right) = \frac{1}{2\sqrt{2}} \left(\begin{array}{c} \sqrt{3} - 1 \\ \sqrt{3} + 1 \end{array} \right)\).

    所以, 可以得出 \(\begin{equation*} \begin{cases} a = d = \frac{1}{2\sqrt{2}}(\sqrt{3} + 1) \\ b = c = \frac{1}{2\sqrt{2}}(\sqrt{3} - 1) \\ \end{cases} \end{equation*}\), 这个变换可以用矩阵 \(M_{r} = \frac{1}{2\sqrt{2}} \left( \begin{array}{c} \sqrt{3} + 1 & \sqrt{3} - 1 \\ \sqrt{3} - 1 & \sqrt{3} + 1 \end{array} \right)\) 表示.

    仅凭这些是没办法看出与代码有什么关系, 这是因为文中代码做了性能优化,

    非常感谢 Kristian Nielsen 的贡献, 不然破头也想不出来.

    先来看一下从 \(xy\) 到 \(uv\) 之间变换: \(\begin{equation*} \begin{cases} u = ax + by \\ v = cx + dy \end{cases} \end{equation*}\), 因为 \(a = d\) 并且 \(b = c\), 所以该关系可以写成 \(\begin{equation*} \begin{cases} u = ax + by \\ v = bx + ay \end{cases} \end{equation*}\).

    整个变换需要 4 次乘法, 文中引入 "Fancy version of zero" 对变换进行重组得到: \(\begin{equation*} \begin{cases} u = ax - bx + bx + by = (a - b)x + b(x + y) \\ v = bx - by + by + ay = (a - b)y + b(x + y) \end{cases} \end{equation*}\),

    其中 \(b(x + y)\) 是重复运算, 可以减少 1 次乘法运算; 再使得 \(a - b = 1\), 则可以再去掉 2 次乘法运算 \(\begin{cases} (a - b)x \\ (a - b)y \end{cases}\),

    那么整个变换只需要 1 次乘法运算 \(b(x + y)\).

    因为 \(a - b = \frac{\sqrt{2}}{2}\), 所以只要让 \(a\) 和 \(b\) 同时除以 \(\frac{\sqrt{2}}{2}\) 或乘以 \(\sqrt{2}\) 即可使得 \(a - b = 1\), 对应了开头说的统一缩放变换,

    可以用矩阵 \(M_{s} = \sqrt{2} \left( \begin{array}{c} 1 & 0 \\ 0 & 1 \end{array} \right)\) 表示.

    同时说明了最终的单形网格空间的基向量并非单位向量: \(\begin{equation*} \begin{cases} \sqrt{2} \times \frac{1}{2\sqrt{2}} = \frac{1}{2} \\ a_{1} = d_{1} = \frac{1}{2}(\sqrt{3} + 1) \\ b_{1} = c_{1} = \frac{1}{2}(\sqrt{3} - 1) \end{cases} \end{equation*}\).

    再回过头来对比整个变换过程 \(M = M_{s}M_{r}\) 和上面的代码, 其中 \(b_{1}(x + y)\) 对应 \(\mathrm{(xin + yin) * F2}\), 并且 \(b_{1} = \mathrm{F2} = \frac{\sqrt{3} - 1}{2}\).

    把单形网格空间上的点还原到正交网格上, 就是求 \(M\) 的逆矩阵 \(M^{-1} = M_{r}^{-1}M_{s}^{-1} = \left( \begin{array}{c} a_2 & b_2 \\ c_2 & d_2 \end{array} \right)\),

    该变换也可以像上面那样写成: \(\begin{equation*} \begin{cases} x = a_{2}u - b_{2}u + b_{2}u + b_{2}v = (a_{2} - b_{2})u + b_{2}(u + v) \\ y = b_{2}u - b_{2}v + b_{2}v + a_{2}v = (a_{2} - b_{2})v + b_{2}(u + v) \end{cases} \end{equation*}\).

    \(M_{s}^{-1} = \frac{\sqrt{2}}{2} \left( \begin{array}{c} 1 & 0 \\ 0 & 1 \end{array} \right)\)

    \(M_{r}^{-1} = \frac{1}{ad - bc} \left( \begin{array}{c} d & -c \\ -b & a \end{array} \right) = \frac{1}{\sqrt{6}} \left( \begin{array}{c} \sqrt{3} + 1 & 1 - \sqrt{3} \\ 1 - \sqrt{3} & \sqrt{3} + 1 \end{array} \right)\)

    \(M^{-1} = M_{r}^{-1}M_{s}^{-1} = \frac{\sqrt{2}}{2} \times \frac{1}{\sqrt{6}} \left( \begin{array}{c} \sqrt{3} + 1 & 1 - \sqrt{3} \\ 1 - \sqrt{3} & \sqrt{3} + 1 \end{array} \right) = \frac{1}{2\sqrt{3}} \left( \begin{array}{c} \sqrt{3} + 1 & 1 - \sqrt{3} \\ 1 - \sqrt{3} & \sqrt{3} + 1 \end{array} \right)\)

    所以, 得出 \(\begin{equation*} \begin{cases} a_{2} = d_{2} = \frac{\sqrt{3} + 1}{2\sqrt{3}} = \frac{3 + \sqrt{3}}{6} \\ b_{2} = c_{2} = \frac{1 - \sqrt{3}}{2 \sqrt{3}} = \frac{\sqrt{3} - 3}{6} \end{cases} \end{equation*}\), 对比代码中的 \(\mathrm{(i + j) * G2}\) 可发现 \(-b_2 = \mathrm{G2} = \frac{3 - \sqrt{3}}{6}\),

    为什么是转负, 因为代码的上下文是这样的:

    double t = (i + j) * G2;
    double X0 = i - t;
    double Y0 = j - t;
    

    \(i\) 和 \(j\) 是单形中第一个被遍历到的角相对坐标 \((0, 0)\) 的分量, \((\mathrm{X0}, \mathrm{Y0})\) 是该单形角对应的正交网格空间坐标, \((\mathrm{x0}, \mathrm{y0})\) 是单形角正交网格空间中上的角指向, 同时是点在单形中的相对坐标.

    到了这里就开始进入下一个环节了: 遍历单形角.

    Ian McEwan 和 Stefan Gustavson 后来又写了 Tiling Simplex Noise and Flow Noise in Two and Three Dimensons, 里面采用了另一种方便计算的二维单形网格, 不过它的单形并非等边的等腰三角形:

    

    这个单形也提高了后续遍历单形角的计算效率.

  • 遍历单形角并求和贡献值
    

    // For the 2D case, the simplex shape is an equilateral triangle.
    // Determine which simplex we are in.
    int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords
    if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1)
    else {i1=0; j1=1;}
    // upper triangle, YX order: (0,0)->(0,1)->(1,1)
    // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
    // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
    // c = (3-sqrt(3))/6
    double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords
    double y1 = y0 - j1 + G2;
    double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords
    double y2 = y0 - 1.0 + 2.0 * G2;
    

    首先是判断点位于哪个单形中, 这里通过点 \(\mathrm{(x0, y0)}\) 在正交网格空间的分量进行判断,

    如果 \(\mathrm{x0} \gt \mathrm{y0}\) 就是菱形的下三角形, 否则就是上三角形.

    因为第二个单形角的相对坐标是 \((0, 1)\) 或 \((1, 0)\), 代码中的 \(\mathrm{x1 = x0 - i1 + G2}\) 完整写法是 \(\mathrm{x0 - [i1 - (1 + 0) \times G2]}\),

    表示第二个单形角在正交网格空间上的角指向 \(x\) 分量, 同理, \(\mathrm{y1}\) 表示对应的 \(y\) 分量.

    因此, \((\mathrm{x1}, \mathrm{y1})\) 是第二个单形角在正交网格空间上的角指向.

    第三个单形角的相对坐标是 \((1, 1)\) 代码中的 \(\mathrm{x2 = x0 - 1 + 2 \times G2}\) 完整写法是 \(\mathrm{x2 = x0 - [1 - (1 + 1) \times G2]}\),

    其中 \(\mathrm{1 - (1 + 1) \times G2}\) 是该单形角的 \(x\) 分量, 所以 \(\mathrm{x2}\) 表示第三个单形角在正交网格空间上的角指向 \(x\) 分量; 同理, \(\mathrm{y2}\) 是对应的 \(y\) 分量.

    因此, \((\mathrm{x2}, \mathrm{y2})\) 是第三个单形角在正交网格空间上的角指向.

    在计算出所有角指向后, 可以开始计算每个角的贡献值了, 下面是原文中的参考代码:

    // Calculate the contribution from the three corners
    double t0 = 0.5 - x0*x0-y0*y0;
    if(t0<0) n0 = 0.0;
    else {
        t0 *= t0;
        n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient
    }
    
    double t1 = 0.5 - x1*x1-y1*y1;
    if(t1<0) n1 = 0.0;
    else {
        t1 *= t1;
        n1 = t1 * t1 * dot(grad3[gi1], x1, y1);
    }
    
    double t2 = 0.5 - x2*x2-y2*y2;
    if(t2<0) n2 = 0.0;
    else {
        t2 *= t2;
        n2 = t2 * t2 * dot(grad3[gi2], x2, y2);
    }
    
    // Add contributions from each corner to get the final noise value.
    // The result is scaled to return values in the interval [-1,1].
    return 70.0 * (n0 + n1 + n2);
    

    这里的 \(\mathrm{n0}\) 正式第一个单形角的贡献值, \(\mathrm{t0}\) 的 4 次方是衰减系数.

    这个衰减系数是用在正交网格空间上的, 正如之前所说的, 衰减函数也代表着一个圆, 它半径为 \(r = \sqrt{0.5}\).

    这个值是如何得到的呢?

    首先在单形网格上求得第一个单形角到它对面边的垂直线, 由于单形是等边三角形, 所以垂直线与对面边的交点正好是对面边的中点, 这样可以确定该边上中点的相对坐标为 \((1, \frac{1}{2})\) 或 \((\frac{1}{2}, 1)\).

    以下三角形为例, 从第一个单形角开始, 往 \(\vec{u}\) 方向增加一个单位, 往 \(\vec{v}\) 方向增加 \(\frac{1}{2}\) 个单位, 最终到达的角相对坐标为 \((1, \frac{1}{2})\).

    如果是上三角形, 最终到达的角相对坐标为 \((\frac{1}{2}, 1)\). 不管哪种情况, 都不影响后续的计算.

    这里以第一个单形角对面边的中点作为例子, 是因为第一个单形角的相对坐标为 \((0, 0)\), 该边中点 \(m_{uv}\) 等于垂直线的向量 \(\vec{m_{uv}} = m_{uv} - (0, 0)\), 垂直线的长度等于 \(|m_{uv}|\).

    这里以下三角形为例, 中点为 \(m_{uv} = (1, \frac{1}{2})\), 接下来把他转换到正交网格空间上, 得到点 \(m_{xy}\):

    \((1 + \frac{1}{2}) \times \mathrm{G2} = \frac{3}{2} \times \frac{3 - \sqrt{3}}{6} = \frac{3 - \sqrt{3}}{4}\)

    \(m_{xy} = (1 - \frac{3 - \sqrt{3}}{4}, \frac{1}{2} - \frac{3 - \sqrt{3}}{4}) = (\frac{1 + \sqrt{3}}{4}, \frac{\sqrt{3} - 1}{4})\)

    最后得到它模长的平方: \(|m_{xy}|^2 = (\frac{1}{4})^2 [(\sqrt{3} + 1)^2 + (\sqrt{3} - 1)^2] = \frac{1}{2}\), 正好就是 \(r^2 = 0.5 = \frac{1}{2}\).

    所以, \(r^2 = 0.5\) 就是单形网格上等边三角形的高在正交网格空间上的长度的平方.

    \(\mathrm{grad3[gi0]}\) 是从预计算好的梯度表格中获取梯度:

    int grad3[][] = {{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},
                     {1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},
                     {0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}};
    

    原文中的二维单形噪声的梯度是三维向量, 而我们在实现时可以完全参考梯度噪声中的梯度计算生成二维梯度向量:

    \(\mathrm{dot(grad3[gi0], x0, y0)}\) 就是梯度的 \(xy\) 分量和角指向之间的点积, \(\mathrm{n1}\) 和 \(\mathrm{n2}\) 的贡献值计算同理.

    private static double dot(int g[], double x, double y) {
        return g[0]*x + g[1]*y;
    }
    

    代码中 \(\mathrm{70 \times (n0 + n1 + n2)}\) 就是单形噪声值, 正如代码的注释所示, \(70\) 是用来把结果归一化成 \([-1, 1]\) 的.

    为什么是 \(70\), 这需要讨论 \(\mathrm{n0 + n1 + n2}\) 的最大值是多少.

    \(\sum\limits_{i=0}^{2} (0.5 - |\vec{d_{i}}|^2)^4 \vec{g_{i}}\vec{d_{i}} = \sum\limits_{i=0}^{2} (0.5 - |\vec{d_{i}}|^2)^4 |\vec{g_{i}}||\vec{d_{i}}|\cos\theta\)

    想让 \(\vec{g_{i}} \cdot \vec{d_{i}}\) 最大, 要满足两个条件:

    1. \(\vec{g_{i}}\) 和 \(\vec{d_{i}}\) 之间的夹角 \(\theta\) 为 0, 换而言之满足这样的关系: \(\vec{g_{i}} = \lambda \vec{d_{i}}\), 其中 \(\lambda \gt 0\);
    2. 两者的模长达到各自的最大值. 由于梯度向量 \(\vec{g}\) 是预设的, 所以可以确定它的最大模长为 \(\sqrt{2}\).

    假设同时满足这两个条件, 上面的等式可以写成 \(\sqrt{2} \sum\limits_{i=0}^{2} (0.5 - |\vec{d_{i}}|^2)^4 |\vec{d_{i}}|\).

    分别讨论三种情况: 输入点位于任意单形角, 位于等边三角形中心以及位于任意边中点.

    这个问题没法直接用费马点(Fermat point)来解决.

    所谓费马点是位于三角形内的一个点, 这个点到三个顶点的距离之和比从其它点算起都要小, 每个三角形只有一个费马点.

    假设三角形的三个顶点分别为 \(A, B, C\), 费马点为 \(P\), 那么距离之和为 \(|\vec{PA}| + |\vec{PB}| + |\vec{PC}|\).

    因为 \(f(x) = (0.5 - x^2)^4 x\) 不是线性关系, 所以 \(\sqrt{2}[f(|\vec{PA}|) + f(|\vec{PB}|) + f(|\vec{PC}|)] \ne \sqrt{2}f(|\vec{PA}| + |\vec{PB}| + |\vec{PC}|)\).

    

    从图中可得知:

    当点处于其中一个角上, 点与对应角的距离为 0, 其贡献值为 \(\sqrt{2} (0.5 - 0)^4 \times 0 = 0\) 其它两个角的贡献值也同样为 0,

    所以 \(\mathrm{n0 + n1 + n2} = 0\);

    当点处于等边三角形的中心, 点到每个角的距离都相等, 中点正好是其高的三等分点, 点到角的距离为 \(\frac{2}{3} \sqrt{\frac{1}{2}} = \frac{\sqrt{2}}{3}\),

    所以 \(\mathrm{n0 + n1 + n2 = 3 \times \sqrt{2}(0.5 - (\frac{\sqrt{2}}{3})^2)^4 \frac{\sqrt{2}}{3}} = 2 (\frac{5}{18})^4 = \frac{625}{52488} \approx \frac{1}{84}\);

    当点处于任意边的中点, 点距离每个角的距离分别为 \(\sqrt{\frac{1}{2}}\), \(\sqrt{\frac{1}{6}}\) 以及 \(\sqrt{\frac{1}{6}}\),

    所以 \(\mathrm{n0 + n1 + n2} = \sqrt{2} [(0.5 - \frac{1}{2})^4 \sqrt{\frac{1}{2}} + (0.5 - \frac{1}{6})^4 \sqrt{\frac{1}{6}} + (0.5 - \frac{1}{6})^4 \sqrt{\frac{1}{6}}] = 2 \sqrt{2} \times \frac{1}{3}^4 \sqrt{\frac{1}{6}} = \frac{2}{81}\sqrt{\frac{1}{3}} \approx \frac{1}{70}\).

    这三个距离中, \(\sqrt{\frac{1}{2}}\) 已经在求 \(m_{xy}\) 时得到了, 另外两个可以根据等边三角形和三角函数得到,

    又或者通过单形网格空间上 \((1, 0)\) 或 \((0, 1)\) 变换到正交网格空间上再进行求模长得出单形的边长, 再除以 2 得到, 这里就不详细说了.

    在对比之下发现, 点处于等边三角形边的中心时贡献值最大.

    然而, 该归一化系数 \(70\) 并不严谨, 目前没找到关于"点处于等边三角形边的中心时贡献值最大"的数学证明, 该结论是通过实验得到的, 后面会提到这个问题.

• 更高维度的单形噪声
  

  从 2D 单噪声的实现中可以看到里面涉及了大量复杂的推导, 那么对于更高维度的单形噪声是不是也要这样呢?

  可以不用, 因为上面推导结果基本上可以在正单形(regular simplex)性质找到.

  这里只讨论几个东西: \(\mathrm{F2}\) 的高维版本 \(F_N\), \(\mathrm{G2}\) 的高维版本 \(G_N\), 正单形的边长 \(L\), 正单形的高 \(H\), 以及更高维度的归一化系数 \(NF\) 如何得到.

  在 \(N\) 维空间上:

  \(F_N = \frac{\sqrt{N + 1} - 1}{N}\)

  \(G_N = \frac{(N + 1) - \sqrt{N+1}}{N^2 + N} = \frac{1 - \frac{1}{\sqrt{N + 1}}}{N}\);

  \(L = \frac{N}{\sqrt{N^2 + N}} = \sqrt{\frac{N}{N + 1}}\);

  这个公式是从 Constructing a regular n-simplex 里面得到的.

  \(H = \sqrt{\frac{N + 1}{2N}} \cdot L = \sqrt{\frac{N + 1}{2N}} \cdot \sqrt{\frac{N}{N + 1}} = \frac{\sqrt{2}}{2}\);

  \(NF(d, r, m, n) = \frac{1}{2m(r^2 - d^2)^n}\), 其中:

  \(d\) 是角指向的模;

  \(r\) 是衰减函数里的半径, 通常是单形角到对边/面的垂直距离;

  \(m\) 是梯度向量的最大模, 为 \(\sqrt{N}\);

  \(n\) 是衰减函数的幂.

  这里的归一化系数的计算方式是基于这个结论得到的: 当点处于单形中的一条边上的中点时, 贡献值达到最大.

  目前我没有找到该支撑该结论的数学证明, 该结论是Sharpe 通过实验得到的 (在评论区里面)得到的.

  需要注意是, Sharpe 计算方式和这里的有些出入, 这里的考虑上了梯度向量的模长不为 1.

  还有, 不管用 Sharpe 的还是说这里的计算方式, 高维 \(NF\) 的计算结果与 Simplex Noise Demystified 的有出入, 考虑到原文中就有不少错误, 因此这里我决定相信 Sharpe.

  这里附上 \(S_N\) 的 Maxima 实现:

  NF(x, r, m, n) := 1 / ((r^2 - x^2)^n * x * m * 2.0);
  

  另外, 为了防止链接失效, 这里还是备份一下评论的内容:

  The way to work out the normalization factor is to think of what the maximum possible value would be generated if given the right conditions. So this would be by taking some sample within the simplex triangle, while having all gradient vectors pointing at that sample. To generate the maximum value it turns out that this sample is halfway along any one of the triangle edges. ( so this means the equal combination of two surflets )

  In 2D, half the edge length of a simplex triangle is…

  x = ( sqrt( 0.5 )/sqrt( 0.75 ) ) * 0.5

  So our 2D normalization factor becomes…

  2DNF = 1.0 / ( x * ( ( 0.5 – x*x ) ^ 4 ) * 2.0 ) = ~99.2043345….

  In 3D, half the edge length of a simplex tetrahedron is…

  x = sqrt( 0.75 ) * 0.5

  So our 3D normalization factor becomes…

  3DNF = 1.0 / ( x * ( ( 0.5 – x*x ) ^ 3 ) * 2.0 ) = ~37.8372272…..

  Something to note: 2D simplex triangles ( as generated by perlins skew/unskew method ) are equilateral. But the 3D simplex tetrahedrons are NOT. They’re slightly skewed…. So the 3D math I’m using here starts to deviate ever so slightly from the actual results. But because it is so minor and to keep math simplicity+sanity, I’ve chosen to ignore this issue.

  (PS: 归一化系数应该是单形噪声中最不严谨的部分了, 在网上找了好几天的相关解释, 发现每个实现的 \(NF\) 都不一样, 着实心累, 考虑到实现中用上近似了, 就不太纠结准确性了.)

细胞噪声

案例一: 估算点与曲线之间的距离

float f( vec2 x ) {
    float r = length(x);
    float a = atan(x.y, x.x);
    return r - 1.0 + 0.5 * sin(3.0 * a + 2.0 * r * r);
}

vec2 grad( vec2 x ) {
    float r = length(x);
    float a = atan(x.y, x.x);
    vec2 da = vec2(-x.y, x.x) / (r * r);
    /* 原文是: vec2 da = vec2(x.y, -x.x) / (r * r), 估计是使用了旋转 90 度的坐标系.
     */
    return (x / r) + (1.5 * da + 2.0 * x) * cos(3.0 * a + 2.0 * r * r);
}

void main() {
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  uv *= 2.0;
  float v = f(uv);
  vec2 g = grad(uv);
  float de = abs(v) / length(g);
  // float de = v / length(g); v 的符号可以用来判断点在图形内还是图形外, 去掉 abs 可以得到一个实心图形
  float lw = 1.0 / 32.0;        // 线粗, 准确来说是线粗占分辨率的百分比, 如果想确定到像素: float lw = 1.0 / iResolution.x;
  float c = smoothstep(0.0, lw, de);
  gl_FragColor = vec4(vec3(c), 1.0);
}

float f( vec2 x ) {
    return dot(x, x) - 0.25;
}

vec2 grad( vec2 x ) {
  return 2.0 * x;
}

void main() {
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  float v = f(uv);
  vec2 g = grad(uv);
  float de = abs(v) / length(g);
  float lw = 1.0 / 32.0;        // 线粗
  float c = smoothstep(0.0, lw, de);
  gl_FragColor = vec4(vec3(c), 1.0);
}

案例二: 估算点与直线之间的距离

float f( vec2 x ) {
  return 3.0 * x.x + 0.1 - x.y;
}

void main() {
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  float v = f(uv);
  vec2 g = vec2(3.0, -1.0);
  float de = abs(v) / length(g);
  float lw = 1.0 / 32.0;        // 线粗
  float c = smoothstep(0.0, lw, de);
  gl_FragColor = vec4(vec3(c), 1.0);
}

初识 SDF

float f2D (vec2 x) { /* ... */ }             // 解析式原函数

vec2 grad2D (vec2 x) { /* ... */ }           // 解析式的导数

float sdf2D (vec2 p) {                       // SDF
  float v = f2D(x);
  float g = grad2D(x);
  return abs(v) / length(g);
}

float f3D (vec3 x) { /* ... */ }             // 解析式原函数

vec3 grad3D (vec3 x) { /* ... */ }           // 解析式的导数

float sdf3D (vec3 p) {                       // SDF
  float v = f3D(x);
  float g = grad3D(x);
  return abs(v) / length(g);
}

对 SDF 描述的图形进行变换

void main() {
  // ...
  vec2 p = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  mat3 trans = mat3(1.0, 0.0, 0.0,
                    0.0, 1.0, 0.0,
                    -0.3, -0.3, 1.0);
  float v = trans * vec3(vec2(p), 1.0);
  float d = sdf2D(v.xy);
  // ...
}

图形组合 - 布尔操作符 (boolean operators)

• 并集 (Union)
  
• 交集 (Intersection)
  
• 差集 (Subtraction)
  
• 异或 (XOR)
  
• 平滑并集 (Smooth Union)
  
• 平滑交集 (Smooth Intersection)
  
• 平滑差集 (Smooth Subtraction)
  

计算点与线段之间的距离

float sdfSegment(vec2 p, vec2 a, vec2 b) {
  vec2 ap = p - a;
  vec2 ab = b - a;
  float ratio = clamp(dot(ap, ab) / dot(ab, ab), 0.0, 1.0);
  return length(ap - ab * ratio);
}

void main() {
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  float c = sdfSegment(uv, vec2(0.3, -0.1), vec2(-0.4, 0.2));
  float l = smoothstep(0.0, 1.0 / iResolution.x, c); // 1px 线粗
  gl_FragColor = vec4(vec3(l), 1.0);
}

计算点与矩形之间的距离

float sdfBox(vec2 p, vec2 size) {
  vec2 d = abs(p) - size;
  return length(max(d, 0.0) + min(max(d.x, d.y), 0.0));
}

void main() {
  #define PI 3.14
  #define ANGLE (PI / 3.0)
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
  // 这里对图形旋转 60 度
  mat2 rotation = mat2(cos(ANGLE), -sin(ANGLE),
                       sin(ANGLE), cos(ANGLE));
  uv = rotation * uv;
  float c = sdfBox(uv, vec2(0.3, 0.1));
  float l = smoothstep(0.0, 1.0 / iResolution.x, c);
  gl_FragColor = vec4(vec3(l), 1.0);
}

计算点与不规则多边形之间的距离

#define N 5

float sdfPoly(vec2 p, vec2[N] vs) {
  const int num = vs.length();
  // 初始化最短距离的平方
  float d = dot(p - vs[0], p - vs[0]);
  float s = 1.0;
  for (int i = 0, j = num - 1; i < num; j=i, i++) {
    vec2 ab = vs[j] - vs[i];
    vec2 ap = p - vs[i];
    float ratio = clamp(dot(ap, ab) / dot(ab, ab), 0.0, 1.0);
    vec2 ao = ap - ab * ratio;

    // 记录最短的距离的平方
    d = min(d, dot(ao, ao));

    bvec3 cond = bvec3( p.y >= vs[i].y, p.y < vs[j].y, ab.x * ap.y > ab.y * ap.x );
    /*
      当 min(vs[i].y, vs[j].y) <= p.y < max(vs[i].y, vs[j].y), 那么射线才和线段 vs[i]vs[j] 相交.

      但是原文并没有这么判断相交, 这是因为还结合了其它条件来判断.

      如果边是往上走的, 那么射线想与边相交, 那么必须满足 vs[i].y < vs[j].y 以及 p.y >= vs[i].y 和 p.y < vs[j].y 同时为 true,

      如果还想要对环绕数进行变化, 那么这时候只能让 cross(ab, ap) > 0, 因为边往上走;

      如果边是往下走, 还想对环绕数进行变化, 那么所有判断必须为 false.
    */
    if ( all(cond) || all(not(cond)) ) s *= -1.0;
  }
  return s * sqrt(d);
}

void main () {
  vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);

  vec2 v0 = 0.8*cos( 0.40 * iTime + vec2(0.0,2.00) + 0.0 );
  vec2 v1 = 0.8*cos( 0.45 * iTime + vec2(0.0,1.50) + 1.0 );
  vec2 v2 = 0.8*cos( 0.50 * iTime + vec2(0.0,3.00) + 2.0 );
  vec2 v3 = 0.8*cos( 0.55 * iTime + vec2(0.0,2.00) + 4.0 );
  vec2 v4 = 0.8*cos( 0.60 * iTime + vec2(0.0,1.00) + 5.0 );

  vec2[] polygon = vec2[](v0, v1, v2, v3, v4);

  float d = sdfPoly(uv, polygon); // [-0.5, 0.5]

  // 分别多边形内外设定颜色
  vec3 col = (d > 0.0) ? vec3(0.9,0.6,0.3) : vec3(0.65,0.85,1.0);

  // 根据距离计算出渐变色的起点色
  col *= 1.0 - exp(-6.0 * abs(d)); // [0.0, 1.0)
  col *= 0.8 + 0.2 * cos(140.0 * d);

  // 进行颜色渐变
  col = mix(col, vec3(1.0), 1.0 - smoothstep(0.0, 0.015, abs(d)));

  gl_FragColor = vec4(vec3(col), 1.0);

  // float l = smoothstep(0.0, 1.0 / iResolution.x, abs(d));
  // gl_FragColor = vec4(vec3(l), 1.0);
}

Ray Marching

• 基础 RayMarching: 计算相机到场景的符号距离场
  

  #define MAX_STEPS 100
  #define MAX_DIST 100.0
  #define SURF_DIST 0.01
  
  // 兼容 glslViewer 和 ShaderToy
  #ifndef PLATFORM_WEBGL
  uniform vec2 u_resolution;
  #define iResolution u_resolution
  #endif
  
  float sphereRadius = 1.0;
  vec3 sphereCenter = vec3(0.0, 1.0, 6.0);
  
  float sdfScene(vec3 p) {
    float sphereDist = length(p - sphereCenter) - sphereRadius;
    float groundDist = p.y;
  
    float d = min(sphereDist, groundDist);
    return d;
  }
  
  float RayMarch(vec3 p0, vec3 rd) {
    float dist = 0.0;
  
    for (int i = 0; i < MAX_STEPS; i++) {
      vec3 p = p0 + dist * rd;
      float d = sdfScene(p);
      dist += d;
      if (dist > MAX_DIST || d < SURF_DIST) break;
    }
  
    return dist;
  }
  
  void main() {
    vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
    vec3 cameraPos = vec3(0.0, 1.0, 0.0);
    vec3 rd = normalize(vec3(uv, 1.0));
    float d = RayMarch(cameraPos, rd);
    d /= 6.0;                     // 避免 d 过大导致画面全白
    gl_FragColor = vec4(vec3(d), 1.0);
  }
  

  

  Figure 29: 3D SDF 效果

  值得注意的是, 整个场景都是通过视点坐标去描述的, 这里的视点坐标系为左手坐标系:

  \(x\) 轴往右为正, \(y\) 轴往上为正, \(z\) 轴往里为正.

  选左手坐标系是因为光线步进渲染啊相比完整的光栅渲染少了投影和 NDC 两个阶段,

  而 NDC 使用的是左手坐标系, 因此作为位于最后一个阶段的视点变换也使用左手坐标系.

  可以通过改变 sphereCenter 以及 cameraPos 进行观察来得到该结论.

• 简单光线
  

  在实现阴影前需要先实现一个简单光线效果: Phong lighting model 的 diffuse 光计算, 具体参考我的另一篇笔记.

  #define MAX_STEPS 100
  #define MAX_DIST 100.0
  #define SURF_DIST 0.01
  
  #ifndef PLATFORM_WEBGL
  uniform vec2 u_resolution;
  uniform float u_time;
  
  #define iResolution u_resolution
  #define iTime u_time
  #endif
  
  float sphereRadius = 1.0;
  vec3 sphereCenter = vec3(0.0, 1.0, 6.0);
  
  float sdfScene(vec3 p) {
    float sphereDist = length(p - sphereCenter) - sphereRadius;
    float groundDist = p.y;
  
    float d = min(sphereDist, groundDist);
    return d;
  }
  
  float RayMarch(vec3 p0, vec3 rd) {
    float dist = 0.0;
  
    for (int i = 0; i < MAX_STEPS; i++) {
      vec3 p = p0 + dist * rd;
      float d = sdfScene(p);
      dist += d;
      if (dist > MAX_DIST || d < SURF_DIST) break;
    }
  
    return dist;
  }
  
  // 新增函数 1
  vec3 GetNormal(vec3 p) {
    // 使用中心差商法估算 ScendSDF 在 p 点上的法线向量
    // 其他计算法线的方法: https://iquilezles.org/articles/normalsSDF/
    vec2 e = vec2(0.01, 0.0);
    float d = sdfScene(p);
    vec3 n = d - vec3(sdfScene(p - e.xyy),
                      sdfScene(p - e.yxy),
                      sdfScene(p - e.yyx));
    return normalize(n);
  }
  
  // 新增函数 2
  float GetLight(vec3 p) {
    // Phong lighting model - diffuse
    mat3 rotateY = mat3(cos(iTime), 0.0, sin(iTime),
                        0.0, 1.0, 0.0,
                        -sin(iTime), 0.0, cos(iTime));
    vec3 lightPos = rotateY * vec3(0.0, 5.0, 6.0); // 对光源进行旋转以便观察效果
    vec3 l = normalize(lightPos - p);
    vec3 n = GetNormal(p);
    float diffuse = clamp(dot(l, n), 0.0, 1.0);
    return diffuse;
  }
  
  void main() {
    vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
    vec3 cameraPos = vec3(0.0, 1.0, 0.0);
    vec3 rd = normalize(vec3(uv, 1.0));
    float d = RayMarch(cameraPos, rd);
  
    // 计算在交点上受到的光照射
    vec3 p = cameraPos + d * rd;
    float diffuse = GetLight(p);
  
    d /= 6.0;                     // 避免 d 过大导致画面全白
    gl_FragColor = vec4(vec3(diffuse), 1.0);
  }
  

  

  Figure 30: 3D SDF 光照

• 阴影实现
  

  有了光线位置就可以计算阴影了,

  #define MAX_STEPS 100
  #define MAX_DIST 100.0
  #define SURF_DIST 0.01
  
  #ifndef PLATFORM_WEBGL
  uniform vec2 u_resolution;
  uniform float u_time;
  
  #define iResolution u_resolution
  #define iTime u_time
  #endif
  
  float sphereRadius = 1.0;
  vec3 sphereCenter = vec3(0.0, 1.0, 6.0);
  
  float sdfScene(vec3 p) {
    float sphereDist = length(p - sphereCenter) - sphereRadius;
    float groundDist = p.y;
  
    float d = min(sphereDist, groundDist);
    return d;
  }
  
  float RayMarch(vec3 p0, vec3 rd) {
    float dist = 0.0;
  
    for (int i = 0; i < MAX_STEPS; i++) {
      vec3 p = p0 + dist * rd;
      float d = sdfScene(p);
      dist += d;
      if (dist > MAX_DIST || d < SURF_DIST) break;
    }
  
    return dist;
  }
  
  vec3 GetNormal(vec3 p) {
    // 使用中心差商法估算 ScendSDF 在 p 点上的法线向量
    vec2 e = vec2(0.01, 0.0);
    float d = sdfScene(p);
    vec3 n = d - vec3(sdfScene(p - e.xyy),
                      sdfScene(p - e.yxy),
                      sdfScene(p - e.yyx));
    return normalize(n);
  }
  
  float GetLight(vec3 p) {
    // Phong lighting model - diffuse
    mat3 rotateY = mat3(cos(iTime), 0.0, sin(iTime),
                        0.0, 1.0, 0.0,
                        -sin(iTime), 0.0, cos(iTime));
    vec3 lightPos = rotateY * vec3(0.0, 5.0, 6.0); // 对光源进行旋转以便观察效果
    vec3 l = normalize(lightPos - p);
    vec3 n = GetNormal(p);
    float diffuse = clamp(dot(l, n), 0.0, 1.0);
  
    /*
      从交点 p 上向光源发射一条阴影射线: float d = RayMarch(p, l);
  
      本应该是这样的, 但根据 RayMarch 函数的实现, p 本来就是交点了, 所以循环的第一轮就因为 d < SURF_DIST 被终止了,
  
      所以我们需要一个相对 p 点稍微偏移一点的点避免这种情况, 只要偏移的地方不是进入物体内, 那么往任意方向偏都可以
    */
    //
    float d = RayMarch(p + n * SURF_DIST * 2.0, l);
    if (d < length(p - lightPos)) {
      diffuse *= 0.1;
    }
  
    return diffuse;
  }
  
  void main() {
    vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
    vec3 cameraPos = vec3(0.0, 1.0, 0.0);
    vec3 rd = normalize(vec3(uv, 1.0));
    float d = RayMarch(cameraPos, rd);
  
    // 计算在交点上受到的光照辐射
    vec3 p = cameraPos + d * rd;
    float diffuse = GetLight(p);
  
    d /= 6.0;                     // 避免 d 过大导致画面全白
    gl_FragColor = vec4(vec3(diffuse), 1.0);
  }
  

  

  Figure 31: 3D SDF 阴影

• 动态计算射线方向以及伽马校正
  

  目前相机的位置以及观察方向都是固定的, 观察方向是 \((0, 0, 1)\).

  并且是在视点空间上的, 可以通过修改 cameraPos 来得到这个结论.

  为了实现相机的控制, 需要先定义函数来构建变换, 实时计算得出视觉射线.

  构建方式类似于 从世界坐标系到视点坐标系的变换 里面矩阵 \(M_{view}\) 的构建过程.

  差别在于, 上述参考是从世界空间到视点空间的变换构建出来的, 这里是直接从视点空间进行构建.

  // 原代码来自: https://www.shadertoy.com/view/3ssGWj
  // 这里对原代码进行了一些调整来提高可读性
  vec3 GetRayDirection(vec2 uv, vec3 cameraPos, vec3 lookAtPos, float near) {
    /* 构建视点坐标系, 这里和上面的参考有一些区别,
  
       在于 GetRayDirection 的参数都是关于左手视点坐标系的
    */
    vec3 up = vec3(0, 1, 0),
      forward = lookAtPos - cameraPos,
      f = normalize(forward),
      right = cross(up, forward),
      r = normalize(right),
      u = cross(f, r);
    /* 因为 f 和 r 都是单位向量且相互垂直,
       u 可以看作是 f 围绕 r 或 r 围绕 f 旋转 90 度到达的方向,
       所以 u 无需再 normalize */
    vec3 nearPlaneCenter = cameraPos + f * near; // 近裁剪平面的中心位置
    // uv 是该平面上的 NDC 坐标, 需要转换到视觉坐标系空间上
    vec3 p = nearPlaneCenter + uv.x * r + uv.y * u;
    vec3 rd = normalize(p - cameraPos);
    return rd;
  }
  
  void main() {
    vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
    vec3 cameraPos = vec3(0.0, 1.0, 0.0);
    vec3 rd = GetRayDirection(uv, cameraPos, cameraPos + vec3(0.0, 0.0, 1.0), 0.5);
    float d = RayMarch(cameraPos, rd);
    vec3 col = vec3(0.0);
  
    if (d < MAX_DIST) {           // 有可能会超出 MAX_DIST
      vec3 p = cameraPos + d * rd;
      float diffuse = GetLight(p);
      col = vec3(diffuse);
    }
  
    col = pow(col, vec3(.4545)); // Gamma correction: [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

  参数 \(near\) 对应着镜头的焦距, 焦距越长, \(fov\) 越小; 反之, 焦距越短, \(fov\) 越大.

  换而言之, \(near\) 可以间接控制 \(fov\).

  另外, 在最后的颜色输出阶段进行了调整, 这项调整被称为伽马校正(gamma correction). 伽马校正是最容易被人忽视的环节.

  Shader 程序一般以线性 RGB (linear RGB)编码像素, 即像素值代表作物理光强度, 且成线性正比关系.

  但像素值在显示器是被解析为电压值, 并且电压值和物理光强度成非线性关系, 比如 2 倍电压不代表 2 倍物理光强度.

  电压值和物理光强度之间为幂函数关系: \(y = x^{\gamma_{\text{monitor}}}\), 称为伽马曲线(gamma curve),

  每个显示器的 \(\gamma_{\text{monitor}}\) 都不一样, 导致它们曲线也不一样, 通常以 \(\gamma_{\text{monitor}} = 2.2\) 作为基准.

  纹理和图片通常以 \((x_{\text{linear}})^{2.2}\) 方式储存像素, 其中 \(x_{\text{linear}}\) 是线性 RGB 编码的像素,

  这种编码方案被称为标准 RGB (sRGB/standard RGB), 所以图片一般是什么颜色, 它在屏幕上显示的就是什么颜色.

  但对于 Shader 来说就不一样了, 把 \(x_{\text{linear}}\) 直接输出到屏幕上, 最终色彩 \((x_{\text{linear}})^{\gamma_{\text{monitor}}}\) 就会偏暗.

  为了输出原本颜色 \(x_{\text{linear}}\), 在被输出到显示器前要根据 \(\gamma_{\text{monitor}}\) 进行校正: \(x^{\gamma_{\text{correction}}}\), 其中 \(\gamma_{\text{correction}} = \frac{1}{\gamma_{\text{monitor}}}\).

  这样在像素输出到浏览器时的色彩为: \([(x_{\text{linear}})^{\gamma_{\text{correction}}}]^{\gamma_{\text{monitor}}} = x_{\text{linear}}\).

  这个处理过程被称为伽马校正(gamma correction).

  

  代码中的 .4545 就是 \(\frac{1}{2.2}\).

  有些渲染引擎会在输出画面前在内部对像素进行一次标准 RGB 编码, 因为它们默认开发者操作的是线性 RGB,

  但这容易出现一个问题: 如果输出贴图的像素作为画面, 那么显示画面就会变暗, 这是因为贴图通常以标准 RGB 编码.

  这种情况需要注意进行伽马校正.

• 使用鼠标控制相机
  

  void main() {
    vec2 uv = (2.0 * gl_FragCoord.xy - iResolution.xy) / min(iResolution.x, iResolution.y);
    vec3 cameraPos = vec3(0.0, 4.0, -5.0); // 重新设置了相机位置
  
    vec2 m = iMouse.xy / iResolution.xy; // iMouse 所处的坐标系也是左下角为原点, 往右边为 x+, 往上为 y+
  
    // 鼠标上移时, 相机往上抬
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    // 鼠标左移时, 镜头左转
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    /* 鼠标动作所对应的变换不一定要和这里一样,
  
       你可以鼠标下移, 相机往上抬, 只需要把这么定义即可:
  
       float sX = sin(-m.y + 0.5);
       float cX = cos(-m.y + 0.5);
  
       你可以鼠标左移, 相机往右转, 只需要把这么定义即可:
  
       float sY = sin(-m.x * 6.2831);
       float cY = cos(-m.x * 6.2831);
     */
  
    cameraPos = rotateY * rotateX * cameraPos;
    // 相机变换是物体变换
  
    vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 3.0, 0.0), 2.0); // 重新设置了 lookAtPos, 以及 near
    float d = RayMarch(cameraPos, rd);
    vec3 col = vec3(0.0);
  
    if (d < MAX_DIST) {           // 有可能会超出 MAX_DIST
      vec3 p = cameraPos + d * rd;
      float diffuse = GetLight(p);
      col = vec3(diffuse);
    }
  
    col = pow(col, vec3(.4545)); // Gamma correction: [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

• 代码结构优化
  

  一旦场景复杂起来, 描述场景的函数 sdfScene 就会难以阅读, 为了让 sdfScene 的结构清晰, 可把物体给封装起来,

  // 球体 SDF
  float sdfSphere (vec3 p, float radius) {
    return length(p) - radius;
  }
  
  // 地面 SDF
  float sdfGround (vec3 p) {
    return p.y;
  }
  
  // sdfScene 就是获取距离 p 点最短的正交距离
  float sdfScene(vec3 p) {
    float sphereDist = sdfSphere(p - sphereCenter, sphereRadius);
    float groundDist = sdfGround(p);
  
    float d = min(sphereDist, groundDist);
    return d;
  }
  

• 抗锯齿 (Anti-Aliasing)
  

  锯齿的本质就是像素之间的变化缺乏连续性, 把连续对象转换成离散对象是造成连续性缺失的原因.

  比如一个像素的颜色为红色, 相邻像素的颜色为蓝色, 两者之间没有一个过渡, 导致看起来很突兀.

  下面使用 盒状模糊 对像素的颜色进行重新计算, 让像素之间的变换更加平滑, 从而实现抗锯齿:

  void main() {
    vec3 cameraPos = vec3(0.0, 4.0, 5.0);
  
    vec2 m = iMouse.xy / iResolution.xy;
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    cameraPos = rotateY * rotateX * cameraPos;
    vec3 col = vec3(0.0);
  
    // AA
    float AASize = 3.0;
    float rAASize = 1.0 / AASize;
    int count = 0;
    for (int aaY = -1; aaY <= 1; aaY++) {
      for (int aaX = -1; aaX <= 1; aaX++) {
        vec2 coord = gl_FragCoord.xy + vec2(aaX, aaY) * rAASize;
        /* 如果 AASize = 1, 那么就相当于把当前像素周围的 8 个像素颜色混合, 使得图像过度模糊,
           所以 AASize 需要适当大小
         */
        vec2 uv = (2.0 * coord - iResolution.xy) / min(iResolution.x, iResolution.y);
        vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 3.0, 0.0), 0.5);
        float d = RayMarch(cameraPos, rd);
        if (d < MAX_DIST) {
          vec3 p = cameraPos + d * rd;
          float diffuse = GetLight(p);
          col += vec3(diffuse);
        }
        count++;
      }
    }
    col /= count;
    // AA ends here
  
    col = pow(col, vec3(.4545)); // Gamma correction: [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

  当然, 抗锯齿会严重影响渲染效率, 效率和效果一直都要进行取舍.

• 无 UV 贴图映射 (Texture Mapping without UVs)
  

  RayMarching 渲染器中不存在顶点数据, 也不存在与顶点关联的纹理坐标.

  即便是光栅渲染器, 也有存在顶点数据当不存在纹理坐标的情况, 需要根据顶点数据计算出纹理坐标.

  如何计算顶点数据取决于具体的几何体以及想要如何映射.

  教材 Fundamentals Of Computer Graphics, 4th, 11.2.1 Geometrically Determined Coordinates 中提供了简单几何体的贴图映射算法.

  最简单贴图映射算法莫过于平面映射(Planar Mapping): 直接让贴图往几何体的一个侧面进行正交投影.

  

  如果从右侧(\(X\) 轴方向)进行投影, 以顶点坐标的 \(yz\) 分量作为 \(uv\);

  如果从顶部(\(Y\) 轴方向)进行投影, 以顶点坐标的 \(xz\) 分量作为 \(uv\);

  如果从正面(\(Z\) 轴方向)进行投影, 以顶点坐标的 \(xy\) 分量作为 \(uv\).

  这种映射算法适用于只展示几何体单面的情况, 几何体的其它面会有严重的拉伸和形变.

  想要解决其它面的拉伸和形变问题, 有很多种方案, 比如很多人第一个想到的就是 3D 贴图,

  3D 贴图也无需 UV 映射, 直接根据顶点坐标进行采样, 不过会占用更高的内存.

  因此, 人们想到了一种比 3D 贴图更低消耗的映射方法: 三平面映射(Triplanar Mapping), 它是平面映射的高阶版本:

  对几何体的顶点 \((x, y, z)\) 分别在 \(XYZ\) 三轴进行正交投影, 分别读取对应的纹理: \(tex_{x}\), \(tex_{y}\) 和 \(tex_{z}\), 以某种权重进行求和, 求和结果就作为该顶点的纹理.

  

  Figure 32: 三平面映射, 图片来源: https://docs.otoy.com

  最简单的权重分配就是求平均值: \((tex_{x} + tex_{y} + tex_{z}) / 3\).

  显然平均值不太准确, 最正式做法应该是根据顶点的法线向量进行权重分配, 正好多元函数的法线向量本质上是梯度向量, 梯度向量由每个坐标轴上的偏导数(变化率)构成.

  具体做法:

  1. 把三个纹理组合成矩阵 \(M = \left( \begin{array}{c} tex_{x} \\ tex_{y} \\ tex_{z} \end{array} \right)\),
  2. 以法线向量 \(n\) 的绝对值 \(abs(n)\) 作为权重, 计算 \(M \cdot abs(n)\), 计算结果作为顶点的纹理.

  以下是具体代码:

  vec3 TriplanarMap(vec3 p, vec3 normal, sampler2D tex) {
    mat3 triMapSamples = mat3(texture(tex, p.yz).rgb,
                              texture(tex, p.xz).rgb,
                              texture(tex, p.xy).rgb);
    return triMapSamples * abs(normal);
  }
  

  接下来会为球体进行贴图, 还需要形成一套为每个不同对象赋予不同贴图的逻辑,

  比如说, 在判断到射线与物体发生相交后, 除了获得相交距离, 还需要知道物体对应哪个贴图.

  为了实现这套逻辑, 整个程序需要进行大幅度的调整, 所以这里贴出整个程序的代码:

  #define MAX_STEPS 100
  #define MAX_DIST 100.0
  #define SURF_DIST 0.01
  // 调整1: 定义空材质ID
  #define MATERIAL_NULL 0
  
  #ifndef PLATFORM_WEBGL
  uniform vec2 u_resolution;
  uniform float u_time;
  uniform vec2 u_mouse;
  uniform sampler2D u_tex0;       // 用于球体的贴图
  
  #define iResolution u_resolution
  #define iTime u_time
  #define iMouse u_mouse
  #define iChannel0 u_tex0
  #endif
  
  // 调整2: 定义结构体以记录交点信息
  struct Intersection {
    float dist;                   // signed distance
    int matID;                    // 材质 ID
  };
  
  float sphereRadius = 1.0;
  vec3 sphereCenter = vec3(0.0, 1.0, 6.0);
  
  float sdfRoundCone( vec3 p, float r1, float r2, float h )
  {
    // sampling independent computations (only depend on shape)
    float b = (r1-r2)/h;
    float a = sqrt(1.0-b*b);
  
    // sampling dependant computations
    vec2 q = vec2( length(p.xz), p.y );
    float k = dot(q,vec2(-b,a));
    if( k<0.0 ) return length(q) - r1;
    if( k>a*h ) return length(q-vec2(0.0,h)) - r2;
    return dot(q, vec2(a,b) ) - r1;
  }
  
  float sdfSphere (vec3 p, float radius) {
    return length(p) - radius;
  }
  
  float sdfGround (vec3 p) {
    return p.y;
  }
  
  // 调整3: 定义并集运算符, 以方便返回 Intersection 类型
  Intersection opU(Intersection a, Intersection b) {
    return a.dist < b.dist ? a: b;
  }
  
  Intersection sdfScene(vec3 p) {
    Intersection sphere = Intersection(sdfSphere(p - sphereCenter, sphereRadius), 1);
    Intersection ground = Intersection(sdfGround(p), 2);
    return opU(sphere, ground);
  }
  
  // 调整4: 修改返回值类型
  Intersection RayMarch(vec3 p0, vec3 rd) {
    Intersection result = Intersection(0.0, MATERIAL_NULL);
  
    for (int i = 0; i < MAX_STEPS; i++) {
      vec3 p = p0 + result.dist * rd;
      // 调整5: 更新交点记录
      Intersection d = sdfScene(p);
      result.dist += d.dist;
      result.matID = d.matID;
      // 调整6: 读取 dist 字段
      if (result.dist > MAX_DIST || d.dist < SURF_DIST) break;
    }
  
    return result;
  }
  
  vec3 GetNormal(vec3 p) {
    vec2 e = vec2(0.01, 0.0);
    // 调整7: 修改返回值类型
    Intersection d = sdfScene(p);
    // 调整8: 读取 dist 字段
    vec3 n = d.dist - vec3(sdfScene(p - e.xyy).dist,
                           sdfScene(p - e.yxy).dist,
                           sdfScene(p - e.yyx).dist);
    return normalize(n);
  }
  
  float GetLight(vec3 p) {
    mat3 rotateY = mat3(cos(iTime), 0.0, sin(iTime),
                        0.0, 1.0, 0.0,
                        -sin(iTime), 0.0, cos(iTime));
    vec3 lightPos = rotateY * vec3(0.0, 5.0, 6.0);
  
    vec3 l = normalize(lightPos - p);
    vec3 n = GetNormal(p);
  
    float diffuse = clamp(dot(l, n), 0.0, 1.0);
  
    // 调整9: 修改返回值类型
    Intersection d = RayMarch(p + n * SURF_DIST * 2.0, l);
    // 调整10: 读取 dist 字段
    if (d.dist < length(p - lightPos)) {
      diffuse *= 0.1;
    }
  
    return diffuse;
  }
  
  vec3 GetRayDirection(vec2 uv, vec3 cameraPos, vec3 lookAtPos, float near) {
    // 构建左手视点坐标系
    vec3 up = vec3(0, 1, 0),
      forward = lookAtPos - cameraPos,
      f = normalize(forward),
      right = cross(up, forward),
      r = normalize(right),
      u = cross(f, r);
    /* 因为 f 和 r 都是单位向量且相互垂直,
       f 旋转到 r 时是围绕 Y+ 以顺时针方向旋转 90 度到达的,
       所以 u 无需再 normalize */
    vec3 nearPlaneCenter = cameraPos + f * near; // 近裁剪平面的中心位置
    // uv 是该平面上的 NDC 坐标, 需要转换到 3D 空间上的坐标
    vec3 p = nearPlaneCenter + uv.x * r + uv.y * u;
    vec3 rd = normalize(p - cameraPos);
    return rd;
  }
  
  vec3 TriplanarMap(vec3 p, vec3 normal, sampler2D tex) {
    mat3 triMapSamples = mat3(texture(tex, p.yz).rgb,
                              texture(tex, p.xz).rgb,
                              texture(tex, p.xy).rgb);
    return triMapSamples * abs(normal);
  }
  
  void main() {
  
    vec3 cameraPos = vec3(0., 1.0, 5.);
  
    vec2 m = iMouse.xy / iResolution.xy;
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    cameraPos = rotateY * rotateX * cameraPos;
    vec3 col = vec3(0.0);
  
    float AA_size = 3.0;
    int count = 0;
    for (int aaY = -1; aaY <= 1; aaY++) {
      for (int aaX = -1; aaX <= 1; aaX++) {
        vec2 coord = gl_FragCoord.xy + vec2(aaX, aaY) / AA_size;
        vec2 uv = (2.0 * coord - iResolution.xy) / min(iResolution.x, iResolution.y);
        vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 0.0, 1.0), 2.0);
        // 调整11: 修改返回值类型
        Intersection d = RayMarch(cameraPos, rd);
        if (d.dist < MAX_DIST) {
          // 调整12: 读取 dist 字段
          vec3 p = cameraPos + d.dist * rd;
          vec3 n = GetNormal(p);
          vec3 tex = vec3(1.0);
          // 调整13: 通过 matID 字段判断射线相交的对象, 为对应对象进行贴图映射
          if (1 == d.matID) {
            tex = TriplanarMap(p, n, iChannel0);
          }
          float diffuse = GetLight(p);
          col += vec3(tex * diffuse);
        }
        count++;
      }
    }
    col /= count;
    col = pow(col, vec3(.4545));  // Gamma correction: [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

  

• 天空盒 (Skybox)
  

  没有与任何物体相交的地方就是天空.

  void main() {
  
    vec3 cameraPos = vec3(0., 1.0, 5.);
  
    vec2 m = iMouse.xy / iResolution.xy;
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    cameraPos = rotateY * rotateX * cameraPos;
    vec3 col = vec3(0.0);
  
    float AA_size = 3.0;
    int count = 0;
    for (int aaY = -1; aaY <= 1; aaY++) {
      for (int aaX = -1; aaX <= 1; aaX++) {
        vec2 coord = gl_FragCoord.xy + vec2(aaX, aaY) / AA_size;
        vec2 uv = (2.0 * coord - iResolution.xy) / min(iResolution.x, iResolution.y);
        vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 0.0, 1.0), 2.0);
        Intersection d = RayMarch(cameraPos, rd);
        if (d.dist < MAX_DIST) {
          vec3 p = cameraPos + d.dist * rd;
          vec3 n = GetNormal(p);
          vec3 tex = vec3(1.0);
          if (1 == d.matID) {
            tex = TriplanarMap(p, n, iChannel0);
          }
          float diffuse = GetLight(p);
          col += vec3(tex * diffuse);
        } else {
          // 天空盒就是没有与射线命中的情况
          col = vec3(0.3, 0.36, 0.6) - (rd.y * 0.4);
        }
        count++;
      }
    }
    col /= count;
    col = pow(col, vec3(.4545)); // Gamma correction: [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

• 雾体 (Fog)
  

  雾体效果的实现方式并不固定, 每个人都可以自行建立雾体的模型, 对于初学者而言, 可以参考 glFog 里面的 GL_EXP 模式的雾体.

  需要注意, 文档里面的公式排版有一些错误, 因此对公式进行了一些调整.

  首先是雾体指数的计算: \(f = e^{-\mathrm{density} \cdot d}\),

  其中 \(d\) 是顶点的视点坐标的 \(z\) 分量, 换句话就是视野距离; \(\mathrm{density}\) 是雾的密度;

  使用 \(C_{f}\) 表示雾体颜色, \(C_{r}\) 表示当前片元颜色, 再根据雾体指数对 \(C_{f}\) 和 \(C_{r}\) 进行混合,

  视野距离越大, 就越接近雾体颜色, \(f\) 是一个随着 \(d\) 增加而减少的函数, 对它"取反": \(F(d) = 1 - f(d)\).

  \(C_{r} \times (1 - F) + F \times C_{f} = C_{r} \times f + (1 - f) \times C_{f}\).

  具体实现如下:

  vec3 applyFog(vec3 fragColor, vec3 fogColor, float d, float start, float density) {
    float fogExp = exp(-(d - start) * density);
    vec3 col = mix(fogColor, fragColor, fogExp);
  
    // 或者
    // float fogExp = 1.0 - exp(-d * density);
    // vec3 col = mix(fragColor, fogColor, fogExp);
    return col;
  }
  
  void main() {
  
    vec3 cameraPos = vec3(0., 4.0, 3.0);
  
    vec2 m = iMouse.xy / iResolution.xy;
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    cameraPos = rotateY * rotateX * cameraPos;
    vec3 col = vec3(0.0);
  
    float AA_size = 3.0;
    int count = 0;
    for (int aaY = -1; aaY <= 1; aaY++) {
      for (int aaX = -1; aaX <= 1; aaX++) {
        vec2 coord = gl_FragCoord.xy + vec2(aaX, aaY) / AA_size;
        vec2 uv = (2.0 * coord - iResolution.xy) / min(iResolution.x, iResolution.y);
        vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 4.0, 0.0), 2.0);
        // 调整1: 把交点计算拿出来
        Intersection d = RayMarch(cameraPos, rd);
        vec3 p = cameraPos + d.dist * rd;
        if (d.dist < MAX_DIST) {
          vec3 n = GetNormal(p);
          vec3 tex = vec3(1.0);
          if (1 == d.matID) {
            tex = TriplanarMap(p, n, iChannel0);
          }
          float diffuse = GetLight(p);
          col += vec3(tex * diffuse);
        } else {
          col = vec3(0.3, 0.36, 0.6) - (rd.y * 0.4);
        }
        // 调整2: 进行雾体计算
        col = applyFog(col, vec3(1., 1., 1.), p.z, 8.0, 1.0 / 80.0);
        count++;
      }
    }
    col /= count;
    col = pow(col, vec3(.4545));  // [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

  例子中球体在原本的位置受雾体的影响不是很大, 可以把球体设置得远一点观察效果.

• 环境光遮蔽 (Ambient Occlusion)
  

  环境光遮蔽用来光线在物体之间相互遮挡产生的阴影效果, 简称 AO.

  在 RayMarching 渲染器中, 通常采用 Alex Evans 的算法来实现 AO.

  为了理解这套算法, 需要了解在 RayMarching 渲染器中如何判断物体自身或物体之间是否存在遮蔽.

  原理很简单, 主要是利用了 RayMarching 渲染器的设计. 观察下图,

  

  其中 \(P\) 点为视觉射线与地面的交点, \(N\) 为地面在 \(P\) 点上的法线向量,

  \(Q\) 是射线从 \(P\) 点出发移动距离 \(d\) 后所到达的点: \(Q = P + d \cdot N\).

  左图只有地面, 右图有地面和物体 \(O\).

  在左图中, \(Q\) 在场景中的距离场 (\(\mathrm{sdfScene}(P + d \cdot N)\)) 远远大于 0;

  在右图中, \(Q\) 在场景中的距离场非常接近 0.

  换而言之, 在 \(Q\) 沿着交点 \(P\) 的法线方向 \(N\) 移动距离 \(d\) 时,

  如果 \(Q\) 的距离场接近于 0 且 \(d\) 非常小, 那么可以断言 \(P\) 被遮蔽了.

  在拥有判断是否遮蔽的方法后, 就得思考如何确定遮蔽的程度.

  再仔细思考可以发现以下规律:

  在没有遮蔽️物的情况下, \(Q\) 的距离场就等于 \(d\): \(d - \mathrm{sdfScene}(P + d \cdot N) = 0\).

  在有遮蔽物的情况下, \(Q\) 的距离场就小于 \(d\): \(d - \mathrm{sdfScene}(P + d \cdot N) \gt 0\).

  正好可以利用 \(d\) 与 \(\mathrm{sdfScene}(P + d \cdot N)\) 的差值作为遮蔽程度.

  在没有遮蔽物的情况下 \(d\) 与 \(\mathrm{sdfScene}(P + d \cdot N )\) 的差值为 0;

  遮蔽程度越高, \(d\) 与 \(\mathrm{sdfScene}(P + d \cdot N)\) 之间的差值越大.

  在具体实现中, 会沿着 \(P\) 点的法线方向以固定步长移动 \(n\) 次, 对每一次的遮蔽程度进行求和, 以此求出 \(P\) 点的遮蔽程度:

  \(\mathrm{occlusion} = \sum \limits_{i=0}^{n}(d_{i} - \mathrm{sdfScene}(P + d_{i} \cdot N))\)

  这里参考 Inigo Quilez 的实现:

  /* 由于参数类型的不同, 这里对原代码进行了一些调整,
     提高了可读性, 但损失了一些性能 */
  float calcAO( in vec3 pos, in vec3 normal, in int sampleCnt )
  {
    float occlusion = 0.0;
    float scale = 1.0;
  
    /* 在实际开发中, 采用数量 sampleCnt 一般是固定的,
       所以会手动计算出 maxOcclusion 然后硬编码 */
    float maxOcclusion = 0.0;
    for (int i=0; i<sampleCnt; i++) {
      maxOcclusion += 0.01 + 0.12 * float(i) / 4.0;
    }
    float iMaxOcclusion = 1.0 / maxOcclusion;
  
    for( int i=0; i<sampleCnt; i++ ) {
      float di = 0.01 + 0.12 * float(i) / 4.0;
      float sdf = sdfScene( pos + di * normal ).dist;
      occlusion += (di - sdf) * scale;
      scale *= 0.95;
      if( occlusion > maxOcclusion ) break;
    }
    // 返回结果是"开放"程度
    return clamp( 1.0 - iMaxOcclusion * occlusion, 0.0, 1.0 );
  }
  
  // 为了展示 AO 效果, 在场景里面添加多一个球体, 使得球体之间的位置接近
  Intersection sdfScene(vec3 p) {
    Intersection sphere1 = Intersection(sdfSphere(p - sphereCenter, sphereRadius), 1);
    Intersection sphere2 = Intersection(sdfSphere(p - (sphereCenter + vec3(0, 2.1, 0)), sphereRadius), 1);
    Intersection ground = Intersection(sdfGround(p), 2);
    Intersection result = opU(sphere1, sphere2);
    result = opU(result, ground);
    return result;
  }
  
  void main() {
  
    vec3 cameraPos = vec3(0., 4.0, 3.0);
  
    vec2 m = iMouse.xy / iResolution.xy;
    float sX = sin(m.y - .5);
    float cX = cos(m.y - .5);
    mat3 rotateX = mat3(1.0, 0.0, 0.0,
                        0.0, cX, sX,
                        0.0, -sX, cX);
  
    float sY = sin(m.x * 6.2831);
    float cY = cos(m.x * 6.2831);
    mat3 rotateY = mat3(cY, 0.0, -sY,
                        0.0, 1.0, 0.0,
                        sY, 0.0, cY);
  
    cameraPos = rotateY * rotateX * cameraPos;
    vec3 col = vec3(0.0);
  
    float AA_size = 3.0;
    int count = 0;
    for (int aaY = -1; aaY <= 1; aaY++) {
      for (int aaX = -1; aaX <= 1; aaX++) {
        vec2 coord = gl_FragCoord.xy + vec2(aaX, aaY) / AA_size;
        vec2 uv = (2.0 * coord - iResolution.xy) / min(iResolution.x, iResolution.y);
        vec3 rd = GetRayDirection(uv, cameraPos, vec3(0.0, 4.0, 0.0), 2.0);
        Intersection d = RayMarch(cameraPos, rd);
        vec3 p = cameraPos + d.dist * rd;
        vec3 n = GetNormal(p);    // 把法线计算拿出来
        float occ = calcAO(p, n, 5); // 计算开放程度, 或称遮蔽修正系数
  
        if (d.dist < MAX_DIST) {
          vec3 tex = vec3(1.0);
          if (1 == d.matID) {
            tex = TriplanarMap(p, n, iChannel0);
          }
          float diffuse = GetLight(p);
          col += vec3(tex * diffuse);
        } else {
          col = vec3(0.3, 0.36, 0.6) - (rd.y * 0.4);
        }
        col = applyFog(col, vec3(1., 1., 1.), p.z, 8.0, 1.0 / 80.0) * occ;  // occ 修正
        count++;
      }
    }
    col /= count;
    col = pow(col, vec3(.4545));  // [0, 1]^3 to [0, 1]^3
    gl_FragColor = vec4(col, 1.0);
  }
  

  

  遮蔽程度的计算并非固定不变的, 对比 Alex Evans 的文章和 Inigo Quilez 的实现后会发现两者在核心原理上保持一致, 但遮蔽计算有差别.

• 包围体层次结构 (Bounding Volume Hierarchy, BVH)
  

  当场景越来越复杂时, 渲染器的渲染效率会大幅度下降, 当然这是无法避免的, 但可以进行一定程度的优化.

  渲染效率下降是由于 sdfScene 的复杂度增加, 比如增加多两个球体:

  Intersection sdfScene(vec3 p) {
    Intersection sphere1 = Intersection(sdfSphere(p - sphereCenter, sphereRadius), 1);
    Intersection sphere2 = Intersection(sdfSphere(p - (sphereCenter + vec3(0, 2.1, 0)), sphereRadius), 1);
  
    // 新增球体
    vec3 transCenter = sphereCenter + vec3(10.0, 0.0, 0.0);
    Intersection sphere3 = Intersection(sdfSphere(p - transCenter, sphereRadius), 1);
    Intersection sphere4 = Intersection(sdfSphere(p - (transCenter + vec3(0, 2.1, 0)), sphereRadius), 1);
  
    Intersection ground = Intersection(sdfGround(p), 2);
  
    Intersection result = opU(sphere1, sphere2);
    result = opU(result, sphere3);
    result = opU(result, sphere4);
    result = opU(result, ground);
    return result;
  }
  

  

  目前场景还不算太复杂, 然而性能也下降不少了.

  在计算其中两个球体的距离场时, 另两个球体的距离场计算是没有必要的, 因为两组球体之间的距离相差比较远, 无意义的计算就是渲染效率下降的原因.

  解决这个问题的有效方法就是按照物体之间的距离把场景划分成若干个区域, 接近的物体会被划分在一个区域,

  区域也可能会以同样方法划分成更多小区域, 整个场景结果最终形成一个树状结构, 比如下图:

  

  这个就是 BVH. 那么问题来了, BVH 是如何提高渲染效率呢?

  其实非常简单, 当射线进入某个区域时, 就只针对这个区域里的物体进行距离场计算, 这样一来就可以避免无意义的距离场计算.

  在实现中可以手动对物体进行划分, 还是以 sdfScene 为例, 可以把球体按照距离进行两两划分.

  // 确定区域的工具
  float sdfBox( vec3 p, vec3 b ) {
    // b = (0.5 * w, 0.5 * h, 0.5 * d)
    vec3 d = abs(p) - b;
    return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0));
  }
  
  Intersection sdfScene(vec3 p) {
    Intersection result = Intersection(sdfGround(p), 2);
    vec3 transCenter = sphereCenter + vec3(10.0, 0.0, 0.0);
  
    /* sdfBox(p - vec3(sphereCenter.x, sphereCenter.y, 2.1), vec3(2.1, 2.1, 4.2)) 是一个盒子,
  
       盒子在位置上与第一个和第二个球体重叠, 在大小上包含了两个球体
    */
  
    if (sdfBox(p - vec3(sphereCenter.x, 2.1, sphereCenter.z), vec3(1.05, 2.1, 1.05)) < result.dist) {
      Intersection sphere1 = Intersection(sdfSphere(p - sphereCenter, sphereRadius), 1);
      Intersection sphere2 = Intersection(sdfSphere(p - (sphereCenter + vec3(0, 2.1, 0)), sphereRadius), 1);
      result = opU(result, sphere1);
      result = opU(result, sphere2);
    }
  
    if (sdfBox(p - vec3(transCenter.x, 2.1, sphereCenter.z), vec3(1.05, 2.1, 1.05)) < result.dist) {
      Intersection sphere3 = Intersection(sdfSphere(p - transCenter, sphereRadius), 2);
      Intersection sphere4 = Intersection(sdfSphere(p - (transCenter + vec3(0, 2.1, 0)), sphereRadius), 2);
      result = opU(result, sphere3);
      result = opU(result, sphere4);
    }
    return result;
  }
  

  这段代码使用了 sdfBox 对物体进行"包裹", 渲染器首先和这些盒子进行相交判断, 在确定盒子相交后再进一步与里面的物体进行相交判断.

  也并非一定要用盒子对物体进行包裹, 也可以使用其它与被包裹的物体形状相近的物体进行包裹.

3D SDF 变换的注意事项

float sdfScene(vec3 p) {
  vec3 sp = vec3(1.0, 3.0, 5.0) * p; // 对球体进行缩小
  float sphereDist = sdfSphere(sp - sphereCenter, sphereRadius); // 距离场 sphereDist 被拉伸了
  float groundDist = sdfGround(p);

  float d = min(sphereDist, groundDist);
  /*
     因为变换对象不是物体本身, 而是射线, 使得距离场被拉伸,

     距离场被拉伸会导致 RayMarch 以过大的步长步进, 从而进入物体内部甚至是穿过物体,

     所以在某些角度观看球体会发现一些问题.
   */
  return d;
}

float sdfScene(vec3 p) {
  vec3 sp = vec3(1.0, 3.0, 5.0) * p; // 对球体进行缩小
  float sphereDist = sdfSphere(sp - sphereCenter, sphereRadius) / 5.0; // 原本缩放的最大分量是 5.0, 所以这里要处以 5.0
  float groundDist = sdfGround(p);

  float d = min(sphereDist, groundDist);
  return d;
}

噪声的导数解析式

// 等同于 mix(a, b, t)
lerp(t, a, b) := a + (b - a) * t;
n(u, v, w) := lerp(w,
  lerp(v, lerp(u, a, b), lerp(u, c, d)),
  lerp(v, lerp(u, e, f), lerp(u, g, h)));
n(u, v, w);