devbox
盐析白兔
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

球谐光照，球谐函数介绍

作者：盐析白兔 | 2024-02-21 21:07:24

赞

踩

球谐函数

///

推导
球谐函数是球面坐标系下的拉普拉斯方程(f(x,y,z) = 0 的泊松方程)中角度部分的解，这里参考 球谐函数 自己手推一遍

从笛卡尔坐标系到球坐标系的推导过程，可以参考：https://wenku.baidu.com/view/f6c75dd2da38376baf1faec3.html

然后主要是分离变量法，考虑

 

虚数域的公式有如下推导，参考 complex analysis - In the real spherical harmonics, where does the sqrt(2) factor come from? - Mathematics Stack Exchange 

比如 UE 里的  SHBasisFunction3，能看到它基本和上面直角坐标系中的表达式一一对应，传进去normal，然后返回出来三阶球谐

FThreeBandSHVector SHBasisFunction3(half3 InputVector)
{
    FThreeBandSHVector Result;
    // These are derived from simplifying SHBasisFunction in C++
    Result.V0.x = 0.282095f; 
    Result.V0.y = -0.488603f * InputVector.y;
    Result.V0.z = 0.488603f * InputVector.z;
    Result.V0.w = -0.488603f * InputVector.x;
 
    half3 VectorSquared = InputVector * InputVector;
    Result.V1.x = 1.092548f * InputVector.x * InputVector.y;
    Result.V1.y = -1.092548f * InputVector.y * InputVector.z;
    Result.V1.z = 0.315392f * (3.0f * VectorSquared.z - 1.0f);
    Result.V1.w = -1.092548f * InputVector.x * InputVector.z;
    Result.V2 = 0.546274f * (VectorSquared.x - VectorSquared.y);
 
    return Result;
}
性质
开始应用之前，首先需要了解一下球谐函数的几个很重要的性质

1.正交完备性

2.旋转不变性 

它的旋转不变不是说怎么转都不变，而是说假定它的旋转变换是 r(x)，我可以先计算旋转再计算球谐 f(r(x))，也可以先计算球谐之后再旋转，也就是  f(r(x)) = r(f(x))

这个的应用一个在于本身采样球是可以旋转的，旋转之后，不用重复计算球谐系数，另外一个在于可以进行空间变换，比如local space 到 world space 的旋转变换，比如有些应用会把球谐系数记录在模型的顶点上，由于是预处理的，一般是local space的，运行时需要转换到 world space，而由于球谐是旋转不变的，因此光照结果不会出错

3.低频高频的问题

经常看到是说球谐适合模拟低频部分的光照信息，比如间接光照，不适合表达高光或者反射等高频的光照信息，对于低频高频，我个人的理解是如下图所示

像图中，所示 m=0 的伴随勒得让多项式 的不同的 l 阶的曲线，阶数越高，频率越高，变化越高，就越能表达更高频的信息(频率高的基底，也只能表达高频信息)，而由于一般实际应用，最多应用3到6阶的球谐，并不足以表达高频信息，所以一般会说球谐适合表达低频信息，高频的可以考虑用小波变换来表示

4.光照函数和传输函数的积分可以转变成它们之间球谐系数的点积，将积分变点乘，大幅简化了计算

传输函数 Tranfer Function 这个概念经常在 PRT(Precompute Radiance Transfer) 提到，我个人的理解是它是对入射照度 Light(x) 进行变换操作的函数，不管是 Lambert 里面的一个 cos，还是标准 BRDF公式 里的 FDG，它们都可以变成传输函数，只是说有些传输函数本身需要的是高频信息，而一般的球谐函数只用有限阶数，适合来表现一些低频信息，不太适合表现高频信息，所以有些函数不会去用球谐表示，而这个特性要做的是说，每种传输函数本来计算都是要组合起来然后积分的，现在都可以各自预计算各自的 球谐系数，然后运行时可以把积分变成系数之间的点积操作，大幅简化了运行时的计算量。

总之有了传输函数，我们就可以表现更加复杂的光照情况，比如遮挡关系等等，而且不会带来很大计算量的开销

应用
了解了基本的性质，就可以开始应用的两步走了，先是投影，然后是重建

先是投影，一般是预计算的，

 

//

另外，这个图不错：Real Spherical Harmonics Figure Table Complex Radial Magnitude - Table of spherical harmonics - Wikipedia

《The Gritty Details》中也列出了0~2阶球谐基函数：

 

 What is the geometric meaning of the inner product of two functions? - Mathematics Stack Exchange

 

至此，你应该理解两条重要结论：

  1. 球谐函数是拉普拉斯方程分离�变量后，角度部分通解的正交项（后面介绍正交性）。
  2. 如何计算球谐函数，可以参见等式（15）（20）（21）。

球谐性质

接着，讨论归一化的球谐函数的性质，它具备两条重要的性质构成了它应用的基石：

  • 正交完备性
  • 旋转不变性

正交完备性

对于任意两个归一化的球谐函数在球面上的积分有

投影的过程就是计算函数积分，计算消耗较大，可以采用离线处理来生成广义傅里叶系数；在实时渲染时，就要简单的线性组合，就可以重建原始函数。当然，由于是有限个系数，就必然存在误差。

我们再来看下连带勒让德函数，如图2所示，随着次数的增加，函数的振动频率会越快。对于函数的展开而言，振动频率越大的基底，它就只能表示越高频的信息，往往一个函数里面的高频信息量是较少的。

 

对于高维矩阵的构造方法非常的复杂，采用的是魏格纳d矩阵（Wigner d-matrices），可以参见文献[4]的讨论，网上也有这个算法的高效实现，有兴趣可以研究研究，参见SHTns。

附录

根据等式（15）的递归关系，就可以很容易计算出连带勒让德函数[3]。

  1. double P(int l,int m,double x)
  2. {
  3.     // evaluate an Associated Legendre Polynomial P(l,m,x) at x
  4.     double pmm = 1.0;
  5.     if(m>0) {
  6.         double somx2 = sqrt((1.0-x)*(1.0+x));
  7.         double fact = 1.0;
  8.         for(int i=1; i<=m; i++) {
  9.             pmm *= (-fact) * somx2;
  10.             fact += 2.0;
  11.         }
  12.     }
  13.     if(l==m) return pmm;
  14.     double pmmp1 = x * (2.0*m+1.0) * pmm;
  15.     if(l==m+1) return pmmp1;
  16.     double pll = 0.0;
  17.     for(int ll=m+2; ll<=l; ++ll) {
  18.         pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
  19.         pmm = pmmp1;
  20.         pmmp1 = pll;
  21.     }
  22.     return pll;
  23. }

根据等式（20）（21），可以计算出球谐函数[3]。

  1. double K(int l, int m)
  2. {
  3.     // renormalisation constant for SH function
  4.     double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
  5.     return sqrt(temp);
  6. }
  7. double SH(int l, int m, double theta, double phi)
  8. {
  9.     // return a point sample of a Spherical Harmonic basis function
  10.     // l is the band, range [0..N]
  11.     // m in the range [-l..l]
  12.     // theta in the range [0..Pi]
  13.     // phi in the range [0..2*Pi]
  14.     const double sqrt2 = sqrt(2.0);
  15.     if(m==0) return K(l,0)*P(l,m,cos(theta));
  16.     else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
  17.     else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
  18. }

在D3D中实现的球谐系数的旋转D3DXSHRotate的实现为：

  1. FLOAT* WINAPI D3DXSHRotate(FLOAT *out, UINT order, const D3DXMATRIX *matrix, const FLOAT *in)
  2. {
  3.     FLOAT alpha, beta, gamma, sinb, temp[36], temp1[36];
  4.  
  5.     TRACE("out %p, order %u, matrix %p, in %p\n", out, order, matrix, in);
  6.  
  7.     out[0] = in[0];
  8.  
  9.     if ((order > D3DXSH_MAXORDER) || (order < D3DXSH_MINORDER))
  10.         return out;
  11.  
  12.     if (order <= 3)
  13.     {
  14.         out[1] = matrix->u.m[1][1] * in[1] - matrix->u.m[2][1] * in[2] + matrix->u.m[0][1] * in[3];
  15.         out[2] = -matrix->u.m[1][2] * in[1] + matrix->u.m[2][2] * in[2] - matrix->u.m[0][2] * in[3];
  16.         out[3] = matrix->u.m[1][0] * in[1] - matrix->u.m[2][0] * in[2] + matrix->u.m[0][0] * in[3];
  17.  
  18.         if (order == 3)
  19.         {
  20.             FLOAT coeff[]={
  21.                 matrix->u.m[1][0] * matrix->u.m[0][0], matrix->u.m[1][1] * matrix->u.m[0][1],
  22.                 matrix->u.m[1][1] * matrix->u.m[2][1], matrix->u.m[1][0] * matrix->u.m[2][0],
  23.                 matrix->u.m[2][0] * matrix->u.m[2][0], matrix->u.m[2][1] * matrix->u.m[2][1],
  24.                 matrix->u.m[0][0] * matrix->u.m[2][0], matrix->u.m[0][1] * matrix->u.m[2][1],
  25.                 matrix->u.m[0][1] * matrix->u.m[0][1], matrix->u.m[1][0] * matrix->u.m[1][0],
  26.                 matrix->u.m[1][1] * matrix->u.m[1][1], matrix->u.m[0][0] * matrix->u.m[0][0], };
  27.  
  28.             out[4] = (matrix->u.m[1][1] * matrix->u.m[0][0] + matrix->u.m[0][1] * matrix->u.m[1][0]) * in[4];
  29.             out[4] -= (matrix->u.m[1][0] * matrix->u.m[2][1] + matrix->u.m[1][1] * matrix->u.m[2][0]) * in[5];
  30.             out[4] += 1.7320508076f * matrix->u.m[2][0] * matrix->u.m[2][1] * in[6];
  31.             out[4] -= (matrix->u.m[0][1] * matrix->u.m[2][0] + matrix->u.m[0][0] * matrix->u.m[2][1]) * in[7];
  32.             out[4] += (matrix->u.m[0][0] * matrix->u.m[0][1] - matrix->u.m[1][0] * matrix->u.m[1][1]) * in[8];
  33.  
  34.             out[5] = (matrix->u.m[1][1] * matrix->u.m[2][2] + matrix->u.m[1][2] * matrix->u.m[2][1]) * in[5];
  35.             out[5] -= (matrix->u.m[1][1] * matrix->u.m[0][2] + matrix->u.m[1][2] * matrix->u.m[0][1]) * in[4];
  36.             out[5] -= 1.7320508076f * matrix->u.m[2][2] * matrix->u.m[2][1] * in[6];
  37.             out[5] += (matrix->u.m[0][2] * matrix->u.m[2][1] + matrix->u.m[0][1] * matrix->u.m[2][2]) * in[7];
  38.             out[5] -= (matrix->u.m[0][1] * matrix->u.m[0][2] - matrix->u.m[1][1] * matrix->u.m[1][2]) * in[8];
  39.  
  40.             out[6] = (matrix->u.m[2][2] * matrix->u.m[2][2] - 0.5f * (coeff[4] + coeff[5])) * in[6];
  41.             out[6] -= (0.5773502692f * (coeff[0] + coeff[1]) - 1.1547005384f * matrix->u.m[1][2] * matrix->u.m[0][2]) * in[4];
  42.             out[6] += (0.5773502692f * (coeff[2] + coeff[3]) - 1.1547005384f * matrix->u.m[1][2] * matrix->u.m[2][2]) * in[5];
  43.             out[6] += (0.5773502692f * (coeff[6] + coeff[7]) - 1.1547005384f * matrix->u.m[0][2] * matrix->u.m[2][2]) * in[7];
  44.             out[6] += (0.2886751347f * (coeff[9] - coeff[8] + coeff[10] - coeff[11]) - 0.5773502692f *
  45.                   (matrix->u.m[1][2] * matrix->u.m[1][2] - matrix->u.m[0][2] * matrix->u.m[0][2])) * in[8];
  46.  
  47.             out[7] = (matrix->u.m[0][0] * matrix->u.m[2][2] + matrix->u.m[0][2] * matrix->u.m[2][0]) * in[7];
  48.             out[7] -= (matrix->u.m[1][0] * matrix->u.m[0][2] + matrix->u.m[1][2] * matrix->u.m[0][0]) * in[4];
  49.             out[7] += (matrix->u.m[1][0] * matrix->u.m[2][2] + matrix->u.m[1][2] * matrix->u.m[2][0]) * in[5];
  50.             out[7] -= 1.7320508076f * matrix->u.m[2][2] * matrix->u.m[2][0] * in[6];
  51.             out[7] -= (matrix->u.m[0][0] * matrix->u.m[0][2] - matrix->u.m[1][0] * matrix->u.m[1][2]) * in[8];
  52.  
  53.             out[8] = 0.5f * (coeff[11] - coeff[8] - coeff[9] + coeff[10]) * in[8];
  54.             out[8] += (coeff[0] - coeff[1]) * in[4];
  55.             out[8] += (coeff[2] - coeff[3]) * in[5];
  56.             out[8] += 0.86602540f * (coeff[4] - coeff[5]) * in[6];
  57.             out[8] += (coeff[7] - coeff[6]) * in[7];
  58.         }
  59.  
  60.         return out;
  61.     }
  62.  
  63.     if (fabsf(matrix->u.m[2][2]) != 1.0f)
  64.     {
  65.         sinb = sqrtf(1.0f - matrix->u.m[2][2] * matrix->u.m[2][2]);
  66.         alpha = atan2f(matrix->u.m[2][1] / sinb, matrix->u.m[2][0] / sinb);
  67.         beta = atan2f(sinb, matrix->u.m[2][2]);
  68.         gamma = atan2f(matrix->u.m[1][2] / sinb, -matrix->u.m[0][2] / sinb);
  69.     }
  70.     else
  71.     {
  72.         alpha = atan2f(matrix->u.m[0][1], matrix->u.m[0][0]);
  73.         beta = 0.0f;
  74.         gamma = 0.0f;
  75.     }
  76.  
  77.     D3DXSHRotateZ(temp, order, gamma, in);
  78.     rotate_X(temp1, order, 1.0f, temp);
  79.     D3DXSHRotateZ(temp, order, beta, temp1);
  80.     rotate_X(temp1, order, -1.0f, temp);
  81.     D3DXSHRotateZ(out, order, alpha, temp1);
  82.  
  83.     return out;
  84. }
  85.  
  86. static void rotate_X(FLOAT *out, UINT order, FLOAT a, FLOAT *in)
  87. {
  88.     out[0] = in[0];
  89.  
  90.     out[1] = a * in[2];
  91.     out[2] = -a * in[1];
  92.     out[3] = in[3];
  93.  
  94.     out[4] = a * in[7];
  95.     out[5] = -in[5];
  96.     out[6] = -0.5f * in[6] - 0.8660253882f * in[8];
  97.     out[7] = -a * in[4];
  98.     out[8] = -0.8660253882f * in[6] + 0.5f * in[8];
  99.     out[9] = -a * 0.7905694842f * in[12] + a * 0.6123724580f * in[14];
  100.  
  101.     out[10] = -in[10];
  102.     out[11] = -a * 0.6123724580f * in[12] - a * 0.7905694842f * in[14];
  103.     out[12] = a * 0.7905694842f * in[9] + a * 0.6123724580f * in[11];
  104.     out[13] = -0.25f * in[13] - 0.9682458639f * in[15];
  105.     out[14] = -a * 0.6123724580f * in[9] + a * 0.7905694842f * in[11];
  106.     out[15] = -0.9682458639f * in[13] + 0.25f * in[15];
  107.     if (order == 4)
  108.         return;
  109.  
  110.     out[16] = -a * 0.9354143739f * in[21] + a * 0.3535533845f * in[23];
  111.     out[17] = -0.75f * in[17] + 0.6614378095f * in[19];
  112.     out[18] = -a * 0.3535533845f * in[21] - a * 0.9354143739f * in[23];
  113.     out[19] = 0.6614378095f * in[17] + 0.75f * in[19];
  114.     out[20] = 0.375f * in[20] + 0.5590170026f * in[22] + 0.7395099998f * in[24];
  115.     out[21] = a * 0.9354143739f * in[16] + a * 0.3535533845f * in[18];
  116.     out[22] = 0.5590170026f * in[20] + 0.5f * in[22] - 0.6614378691f * in[24];
  117.     out[23] = -a * 0.3535533845f * in[16] + a * 0.9354143739f * in[18];
  118.     out[24] = 0.7395099998f * in[20] - 0.6614378691f * in[22] + 0.125f * in[24];
  119.     if (order == 5)
  120.         return;
  121.  
  122.     out[25] = a * 0.7015607357f * in[30] - a * 0.6846531630f * in[32] + a * 0.1976423711f * in[34];
  123.     out[26] = -0.5f * in[26] + 0.8660253882f * in[28];
  124.     out[27] = a * 0.5229125023f * in[30] + a * 0.3061861992f * in[32] - a * 0.7954951525f * in[34];
  125.     out[28] = 0.8660253882f * in[26] + 0.5f * in[28];
  126.     out[29] = a * 0.4841229022f * in[30] + a * 0.6614378691f * in[32] + a * 0.5728219748f * in[34];
  127.     out[30] = -a * 0.7015607357f * in[25] - a * 0.5229125023f * in[27] - a * 0.4841229022f * in[29];
  128.     out[31] = 0.125f * in[31] + 0.4050463140f * in[33] + 0.9057110548f * in[35];
  129.     out[32] = a * 0.6846531630f * in[25] - a * 0.3061861992f * in[27] - a * 0.6614378691f * in[29];
  130.     out[33] = 0.4050463140f * in[31] + 0.8125f * in[33] - 0.4192627370f * in[35];
  131.     out[34] = -a * 0.1976423711f * in[25] + a * 0.7954951525f * in[27] - a * 0.5728219748f * in[29];
  132.     out[35] = 0.9057110548f * in[31] - 0.4192627370f * in[33] + 0.0624999329f * in[35];
  133. }
  134.  
  135. FLOAT * WINAPI D3DXSHRotateZ(FLOAT *out, UINT order, FLOAT angle, const FLOAT *in)
  136. {
  137.     UINT i, sum = 0;
  138.     FLOAT c[5], s[5];
  139.  
  140.     TRACE("out %p, order %u, angle %f, in %p\n", out, order, angle, in);
  141.  
  142.     order = min(max(order, D3DXSH_MINORDER), D3DXSH_MAXORDER);
  143.  
  144.     out[0] = in[0];
  145.  
  146.     for (i = 1; i < order; i++)
  147.     {
  148.         UINT j;
  149.  
  150.         c[i - 1] = cosf(i * angle);
  151.         s[i - 1] = sinf(i * angle);
  152.         sum += i * 2;
  153.  
  154.         out[sum - i] = c[i - 1] * in[sum - i];
  155.         out[sum - i] += s[i - 1] * in[sum + i];
  156.         for (j = i - 1; j > 0; j--)
  157.         {
  158.             out[sum - j] = 0.0f;
  159.             out[sum - j] = c[j - 1] * in[sum - j];
  160.             out[sum - j] += s[j - 1] * in[sum + j];
  161.         }
  162.  
  163.         if (in == out)
  164.             out[sum] = 0.0f;
  165.         else
  166.             out[sum] = in[sum];
  167.  
  168.         for (j = 1; j < i; j++)
  169.         {
  170.             out[sum + j] = 0.0f;
  171.             out[sum + j] = -s[j - 1] * in[sum - j];
  172.             out[sum + j] += c[j - 1] * in[sum + j];
  173.         }
  174.         out[sum + i] = -s[i - 1] * in[sum - i];
  175.         out[sum + i] += c[i - 1] * in[sum + i];
  176.     }
  177.  
  178.     return out;
  179. }

/

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/盐析白兔/article/detail/126567
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号