前言

在之前我们了解了如何使用烘焙贴图的方式计算IBL的diffuse、specular，但这种方式存在些许问题：

  • 纯静态

  • 显存与带宽开销大、ALU

    HDR cubemap 贴图一般是很大的，会占用几十mb

    采样Cubemap带来的开销

本文将细细讲解解决这一问题的球谐函数，它究竟是什么，为什么需要它，该如何实现它，它有什么优缺点

球谐函数的what、why

如果去网上搜索关于球谐函数的讲解，会发现大部分都是往你脸上丢非常复杂的公式，在这就被劝退了。学习，是一个先理解本质与作用，再考虑如何实现的过程，这里不会先看公式，而是说明球谐函数的本质是什么

  • 球谐函数的本质

    球谐函数其实就是在球面上的、由多个sin、cos函数组成的基函数集合（基函数即sin、cos）

    而其中由多个sin、cos函数组成的基函数集合也就是网上所说的傅里叶变换

    那么球谐函数可以理解为，三维空间中二维平面上的傅里叶变换

  • 为什么需要球谐函数

    用于拟合特定曲线（形状）

    为什么球谐函数可以拟合任一曲线呢？这一特性正是sin、cos支持的，只要sin、cos够多，可以拟合任何曲线，初音未来都没问题

    对于IBL的diffuse信息，不难发现是低频信息，如果可以用球谐函数来拟合这些低频信息的形状，那不就可以省去烘焙了嘛。没错，球谐函数正是做这件事的

如何在IBL中运用球谐函数

  • 工业界使用的公式IBL Diffuse如下：

    • 在工业界并不会直接用sin、cos拟合，因为GPU端sin、cos三角函数指令开销是比较高的，更不必说会积分，所以工业界的方案是将三角函数代替为normal的x、y、z，并配合一些炮轰出的magic number
    void EvaluateSHBasis(float3 dir, out float Y[9])
{
    float x = dir.x;
    float y = dir.y;
    float z = dir.z;

    // L0
    Y[0] = 0.282095f;

    // L1
    Y[1] = 0.488603f * y;
    Y[2] = 0.488603f * z;
    Y[3] = 0.488603f * x;

    // L2
    Y[4] = 1.092548f * (x * y);
    Y[5] = 1.092548f * (y * z);
    Y[6] = 0.315392f * (3.0f * z * z - 1.0f);
    Y[7] = 1.092548f * (x * z);
    Y[8] = 0.546274f * (x * x - y * y);
}
    • 本来球面的立体角以公式中的sinθdθdΦ是可以表示的，但实际计算中模型不是球面，而是立方体，因此实际计算中会采用一般曲面积分公式来表示立体角——dΩ = \frac{dA·cosθ}{r^2}，否则会导致某些像素被重复采样了好几次，而某些边角的高频像素可能被彻底漏掉
    // 映射到[-1, 1]
float u_coord = (samplePos.x + 0.5f) * g_SkyboxSize.z * 2.0f - 1.0f;
float v_coord = (samplePos.y + 0.5f) * g_SkyboxSize.w * 2.0f - 1.0f;

// 每次在特定的面计算像素坐标，那么xyz中有一个abs一定为1
// 假设现在是以z = 1平面，采样点坐标为[u, v, 1]，原点到采样点的距离平方r^2 = x^2 + y^2 + 1
float distSq = 1.0f + u_coord * u_coord + v_coord * v_coord;

// θ是V与V夹角,N = (0,0,1)，V = (u, v, 1)，根据点积：N⋅V=∣N∣∣V∣cosθ => cosθ = 1 / r
// 因为归一化了, dA = 2 * 2 = 4
float solidAngle = 4.0f / (sqrt(distSq) * distSq);
    • 为了更好的性能，Lambert 漫反射的公式会除以PI，因此工业界直接对E(n)的系数都除以PI
    finalSH[0] *= 0.28209479f;
finalSH[1] *= 0.32573501f;
finalSH[2] *= 0.32573501f;
finalSH[3] *= 0.32573501f;
finalSH[4] *= 0.27313711f;
finalSH[5] *= 0.27313711f;
finalSH[7] *= 0.27313711f;
finalSH[6] *= 0.07884789f;
finalSH[8] *= 0.13656855f;
    • 在最后我们会加权平均（积分的总立体角可能不会等于4PI），并乘上球体的总立体角表面积4PI

优化

不难发现，如果skybox是动态的，需要每帧或多帧分摊实时计算SH系数，如何优化SH系数的计算过程将非常重要

如果在compte shader里一次性循环skybox的六个面、宽、高，相当于使用一个线程完成了所有计算（所有线程跑相同计算），性能是非常差的

  • 可以想到的方案是将三个循环分摊给所有线程——Dispatch（宽，高，面） 
    m_pCommand->Dispatch(CeilDivide(skyboxTex.GetWidth(), threadGroupSize.x),
                               CeilDivide(skyboxTex.GetHeight(), threadGroupSize.y),
                               6);

但这会有一个问题：每64线程为一组计算SH、weight，需要GroupMemoryBarrierWithGroupSync同步，也就是说需要用到LDS，将每组每个线程得到的SH、weight累加到LDS。但这不是最重要的，最重要的是，LDS的64个数据，如何充分利用到线程？如果是int类型，可以使用InterlockedAdd（），在同步后立即计算，但很不幸这里是float类型，且InterlockedAdd会打破线程的并行，退化为串行，且每个线程都需要跑到全局显存，效率非常低

  • 一种可行的方案：每一个线程组的第0号线程，统一求和64个数据 
    if (localIdx == 0)
{
  float4 totalSH[9] = (float4[9])0;
  float totalWeight = 0.f;

  for (int k = 0; k < 64; ++k)
  {
      for (int j = 0; j < 9; ++j)
      {
          totalSH[j] += g_SHCoefficients[j][k];
      }
      totalWeight += g_SHWeights[k];
  }

  const UINT saveIndex = groupID.x +
                         groupID.y * g_SHCoefficientsTempCount.x +
                         groupID.z * g_SHCoefficientsTempCount.x * g_SHCoefficientsTempCount.y;
  SHCoefficientData groupResult;
  for (int i = 0; i < 9; ++i)
  {
      groupResult.SHCoefficients[i] = totalSH[i];
  }
  groupResult.TotalWeight = totalWeight;
  groupResult._Padding = float3(0, 0, 0);

  Elysia_Save_Temp_SHCoefficient_Data(saveIndex, groupResult);
}

    这种方案不会因为锁等待派对，不会每个线程都去全局显存读取，效率会提升不少

  • 但每组依然会有63个线程无所事事，尤其暴力循环，不如来个二分法，每次折半叠加——如64个线程，每次取一半的线程（32个），前一半的线程叠加后一半的线程，继续折半，直至1，时间复杂度从O(N)降到了O(logN)

    for (UINT stride = 32; stride >= 1; stride >>= 1)
{
  if (localIdx < stride)
  {
      for (UINT i = 0; i < 9; ++i)
      {
          g_SHCoefficients[i][localIdx] += g_SHCoefficients[i][localIdx + stride];
      }
      g_SHWeights[localIdx] += g_SHWeights[localIdx + stride];
      GroupMemoryBarrierWithGroupSync();
  }
}

最终效果如下：

他们曾如此骄傲的活过，贯彻始终