球谐光照——简单理解和实现
转自：https://blog.****.net/qq_33999892/article/details/83862583
1.简介：
游戏中的光照可以被简单分为两个部分，直接光和间接光。
直接光中我们研究各种不同的BRDF材质，甚至BSDF，BSSSRDF等等。这些模型据有很不错的表现力，足够我们区分金属，皮肤，木头等等不同物体的着色表现。
但这并不能满足我们，因为光并不是那么简单，光会被反射，会被折射，会被透射，会被吸收，所以物体的受光情况同时又由这个场景的其他物体决定，这部分光照同时拥有着更加富强的表现力，被我们称作间接光。
用来做这部分光照的算法也有很多，如raytracing，photonmap之流。
这些算法很复杂也很有趣，不过这些并非我们讨论的重点。因为我们进一步思考就会发现，虽然raytracing这类算法已经从场景中的光源出发，一步一步的追踪光线，这已经十分符合我们在物理学光线发生的那样，但是，真实世界环境中的光源是极度复杂的，事实上——你不知道是个什么鬼东西的玩意就在发着微弱的光并且微弱的影响着场景。也就是我们游戏中的光源复杂度完全无法达到现实中的要求。
于是人们提出了一种叫环境光的概念——一个极度复杂的场景被360度采样后保留在一个环境球上。
大概看上去长这样，有时也会直接作为天空盒，用来表现一些游戏的远景。我们寄希望于用他来还原出这个场景会对某个具体的模型产生的光照效果——当然，由于环境光中只保留了方向性，所以这只是间接光计算中的一种辅助手段。
2.解决:
那么，我们如何解决环境光对物体的影响呢？
首先，我们要理解这样一个积分式:
由于环境光我们只考虑方向性，所以在求解一个点的着色时，我们对这个点的方向上做积分：从四面八方来的光(Lenv(s))->在这个方向上没有被遮挡(V(p->s))->由于辐射通量符合余弦衰减(H(s)就是点乘项)->与物体作用发生BRDF的影响(fr(p,s->d))最终积分累加得到的这个点的着色。
2.1球谐光照与预计算：
理解完积分式后，我们很难理解球谐函数在这里的提出，因为我们很自然的会想到，只需要离线先为模型的每个点计算一下这个积分，并存储在文件里，在实时绘制的时候拿出来运用就可以了。
但是这里的困境在于，虽然受到烘焙的物体常常是固定不动的，但是游戏中的天空盒常常我们希望是动态的，我们很自然的会希望在天气系统中，天空盒会动态改变，甚至进行一定的旋转。或许不同的天空盒，我们还可以通过，烘焙多次来解决，但是不同角度的天空盒，总是不能通过每个角度都来烘焙一次，然后实时绘制时用插值的方式来解决。
所以之前的研究中就想到，如果能分离Lenv(s)这一项的积分和其他项就好了。因为这样，与光源相关的积分我们可以单独烘焙，与模型相关的积分我们单独烘焙，到了实时，我们只需加载不同的天空盒与模型的数据组合就好了。
2.1.1球谐基函数:
为了解决这个问题，我们不得不暂时跳出这个积分式，来点数学知识。
在向量中，我们很好理解基向量这一概念 ,其实可以理解为i,j,k是坐标轴，x的坐标为(a,b,c)
在基函数中，我们扩展了这个概念，对于函数，亦将f1,f2,f3称为基函数。我们暂且不讨论基函数其他方面的意义，回归到我们的问题，我们运用他，可以获得怎样好的性质？
也就是对于我们的积分，我们如果将光源函数看做f(s),其余项看做一个统一的k(s)——也被称做预计算传输函数，然后我们将光源函数分解成一组基函数来表达，就可以把ci提出去，这样就做到了光源函数和其余项积分的解耦——ci为光源函数的积分，yi(s)*k(s)为其余项的积分。
特别的说明一下，，此过程又被称为投影。
这里我们有非常多的基函数可以选择，最典型的有球谐函数，小波函数，球高函数。
以球谐函数为例，其形式为:
其中的为球坐标的表达，也就是这个是一个只和方向有关的函数，正好用于我们光源函数的分解。
再来解释一下l,m两项，这里l,m对应着不同的球谐基函数，通常l被称为带，而我们定义第i个球谐基函数为第l*(l+1)+m个。
下面这副图将解释的一清二楚。
可以看到，所有的球谐基函数好像被l划分到不同的带上，而前l带上正好具有l^2个球谐基函数。
2.1.2积分式求解
我们再次来看看原来的积分式——
首先，我们只考虑Diffuse项，所以对于Diffuse项，fr(s,d)也就是BRDF项可以暂时考虑为常数。
积分式变成了——
由2.1.2可知上式等价于
核心问题变成了如何解决两个积分的问题。
这里我们引入蒙特卡洛积分法。
想要看详细推导可以在之后的一些相关链接找到，这里，我们直接使用他的结论，即寻找一个采样的概率函数p(x)，这样就有:
简单来说就是通过把他理解为概率密度函数求期望的过程来求解积分的数值解。
这里我们是对一个球的方向做积分，所以很容易想明白要在球上做均匀的采样，由于球面积是4*pie*r，由于方向与r无关，不妨设r=1，这时可以得到p(x) = 1/4*pie。
于是我们有:
这样就可以很轻松的编程求解积分的数值解。
2.1.3 内反射
上述的积分式我们暂时只考虑了一次反射的问题，但是实际情况中是会发生两次甚至三次反射的，或者说，在游戏里，我们可以根据我们所需来考虑他到底需要几次反射。
那么问题在于，我们如何进行第两次反射的积分？
一个朴素的想法在于——我们可以把第0次反射的点(直接被照亮的点)看作光源。
这个积分式和之前自阴影的差距主要在于，这次是可见性为不可见的时候才能被积。
我们继续推导他——
其中sh0(s)就是预计算传输函数得出的结果，是一组与顶点相关的系数。
接着我们又可以将L1这里的积分预计算出来并和sh0相加得到sh1,这样在实时绘制中，我们不管是需要内反射or自阴影的效果，其实计算规则都是一样的，也就是和ci项的系数相乘后累加即可。
2.1.4 实时绘制与旋转不变性
终于我们完成了所有的预计算，那么实际在实时绘制的时候，其本质只需要把在前一步中被我们解耦的两组数据加载进来，再次相乘并累加即可。但是这里我们不得不提到一个球谐光照的一个有趣性质，旋转不变性。也就是说，他支持你在实时中旋转光源——也就是那个环境球。这是由于这样一组等价：
没错！你可以把原本需要旋转完再进行预计算的事变成一次性SH计算完，根据实时的Rotate角度来对SH系数做变换得到，这样的效果是完全等价的。
也就是我们现在的实时绘制步骤变成了——
1.获得SH光照系数和与模型点相关的传输数据系数
2.计算实时旋转角度
3.旋转SH光照系数
4.旋转后的SH光照系数与传输数据系数相乘累加获得每点的光照颜色值。
这里的难点在于，我们的旋转矩阵通常旋转的是三维齐次坐标，是一个4x4的矩阵，而对SH系数旋转，显然不是用这个矩阵——因为如果你使用的SH系数是16维的话(这是一个游戏中使用的最多的维度)，你需要一个16x16的矩阵才可以。
当然，这个SH旋转矩阵和我们用的4x4旋转矩阵是有一个数学上的对应关系的，但由于他过于复杂，我打算在实现篇直接给出代码。
3.实现
3.1 如何产生采样点？
private Vector2[] MCIntegrator(int maxCount)
{
Vector2[] res = new Vector2[maxCount];
for (int i = 0; i < maxCount; i++)
{
res [i].x = Random.Range (0.0f,1.0f);
res [i].y = Random.Range (0.0f,1.0f);
res [i].x = 2 * Mathf.Acos (Mathf.Sqrt(1-res[i].x));
res [i].y = 2 * res [i].y * Mathf.PI;
}
return res;
}
这里的关键在于，我们先产生一个正方形区域的均匀点，然后将其映射到球面上——这样就可以在球上均匀的采样。其实这个做法本质非常像UV纹理映射的关系——UV是正方形区域，球面是你要映射的纹理Mesh面。（PS：这也是面试的常考题——如何在球上均匀采样？）
3.2 预计算积分的一个大致代码片段：
private void preCompute(int i)
{
Vector3 position = SHMeshObject.vertices[i];
Vector3 normal = SHMeshObject.normals[i];
Vector2[] randomVec = MCIntegrator (maxSamples);
float[] result = new float[maxPower];
for (int j = 0; j < maxPower; j++)
result [j] = 0.0f;
for(int j=0;j<maxSamples;j++)
{
Vector3 tempDir = tranfer (randomVec[j].x,randomVec[j].y);
float csn = Mathf.Clamp01(Vector3.Dot(tempDir.normalized,normal.normalized));
float shadow = RayCast (position,tempDir);
for (int k = 0; k < maxPower; k++)
{
int l = 0;
int m=0;
SHFunc.tranferK (k,ref l,ref m);
float y = SHFunc.SH (l,m,randomVec[j].x,randomVec[j].y);
result [k] += y * shadow * csn*albedo;
}
}
float factor = 4.0f * Mathf.PI / maxSamples;
for (int j = 0; j < maxPower; j++)
{
AllIntegratorTasks [i].coffeeSHResult [j] = result [j] * factor;
}
}
忠实的计算积分即可，一些比较重要的函数如光线投射RayCast,由于这里是在Unity下，直接使用Unity的API——Physics.raycast即可，而如果想使用GL或DX去写，免不了要自己手动写一些类似ocTree的东西对raycast做优化。
3.3SH旋转
public class PermutedMatrix
{
public PermutedMatrix(Matrix4x4 m)
{
mMat44 = m;
}
static int permute(int v)
{
if (v == 1)
return 0;
if (v == -1)
return 1;
if (v == 0)
return 2;
return 0;
}
public float GetByMN(int m,int n)
{
int row = permute (m);
int column = permute (n);
if (row == 0 && column == 0)
return mMat44.m00;
if (row == 0 && column == 1)
return mMat44.m01;
if (row == 0 && column == 2)
return mMat44.m02;
if (row == 1 && column == 0)
return mMat44.m10;
if (row == 1 && column == 1)
return mMat44.m11;
if (row == 1 && column == 2)
return mMat44.m12;
if (row == 2 && column == 0)
return mMat44.m20;
if (row == 2 && column == 1)
return mMat44.m21;
if (row == 2 && column == 2)
return mMat44.m22;
return -1;
}
private Matrix4x4 mMat44;
}
public class SHRotate
{
private static float delta(int m, int n)
{
return (m == n ? 1 : 0);
}
private static void uvw(int l,int m,int n, ref float u,ref float v,ref float w)
{
float d = delta(m, 0);
int abs_m = Mathf.Abs(m);
float denom;
if (Mathf.Abs(n) == l)
denom = (2 * l) * (2 * l - 1);
else
denom = (l + n) * (l - n);
u = Mathf.Sqrt((l + m) * (l - m) / denom);
v = 0.5f * Mathf.Sqrt((1 + d) * (l + abs_m - 1) * (l + abs_m) / denom) * (1 - 2 * d);
w = -0.5f * Mathf.Sqrt((l - abs_m - 1) * (l - abs_m) / denom) * (1 - d);
}
private static float P(int i,int l,int a,int b,PermutedMatrix R,SHRotateMatrix M)
{
if (b == -l)
{
return (R.GetByMN(i,1) * M.GetValueByBand(l - 1, a, -l + 1) + R.GetByMN(i, -1) * M.GetValueByBand(l - 1, a, l - 1));
}
else if (b == l)
{
return (R.GetByMN(i, 1) * M.GetValueByBand(l - 1, a, l - 1) - R.GetByMN(i, -1) * M.GetValueByBand(l - 1, a, -l + 1));
}
else
{
return (R.GetByMN(i, 0) * M.GetValueByBand(l - 1, a, b));
}
}
private static float U(int l,int m,int n,PermutedMatrix R,SHRotateMatrix M)
{
if (m == 0)
return (P(0, l, 0, n, R, M));
return (P(0, l, m, n, R, M));
}
private static float V(int l,int m,int n,PermutedMatrix R,SHRotateMatrix M)
{
if (m == 0)
{
float p0 = P(1, l, 1, n, R, M);
float p1 = P(-1, l, -1, n, R, M);
return (p0 + p1);
}
else if (m > 0)
{
float d = delta(m, 1);
float p0 = P(1, l, m - 1, n, R, M);
float p1 = P(-1, l, -m + 1, n, R, M);
return (p0 * Mathf.Sqrt(1 + d) - p1 * (1 - d));
}
else
{
float d = delta(m, -1);
float p0 = P(1, l, m + 1, n, R, M);
float p1 = P(-1, l, -m - 1, n, R, M);
return (p0 * (1 - d) + p1 * Mathf.Sqrt(1 - d));
}
}
private static float W(int l,int m,int n,PermutedMatrix R,SHRotateMatrix M)
{
if (m == 0)
{
return (0);
}
else if (m > 0)
{
float p0 = P(1, l, m + 1, n, R, M);
float p1 = P(-1, l, -m - 1, n, R, M);
return (p0 + p1);
}
else // m < 0
{
float p0 = P(1, l, m - 1, n, R, M);
float p1 = P(-1, l, -m + 1, n, R, M);
return (p0 - p1);
}
}
private static float M(int l,int m,int n,PermutedMatrix R,SHRotateMatrix M)
{
// First get the scalars
float u=0.0f, v=0.0f, w=0.0f;
uvw(l, m, n, ref u, ref v, ref w);
// Scale by their functions
if (u!=0.0f)
u *= U(l, m, n, R, M);
if (v!=0.0f)
v *= V(l, m, n, R, M);
if (w!=0.0f)
w *= W(l, m, n, R, M);
return (u + v + w);
}
public static Vector3[] Rotate(Vector3[] src,Matrix4x4 rot)
{
SHRotateMatrix shrm = transfer (rot,(int)Mathf.Sqrt(src.Length));
Vector3[] dest = shrm.Transform (src);
return dest;
}
public static SHRotateMatrix transfer(Matrix4x4 rot,int bands)
{
SHRotateMatrix result = new SHRotateMatrix (bands*bands);
result.SetValue (0, 0, 1);
PermutedMatrix pm = new PermutedMatrix(rot);
for (int m = -1; m <= 1; m++)
for (int n = -1; n <= 1; n++)
result.SetValueByBand (1,m,n,pm.GetByMN(m,n));
for (int band = 2; band < bands; band++)
{
for (int m = -band; m <= band; m++)
for (int n = -band; n <= band; n++)
result.SetValueByBand(band,m,n,M(band, m, n, pm,result));
}
return result;
}
}
public class SHRotateMatrix
{
public Vector3[] Transform(Vector3[] src)
{
int bands = (int)Mathf.Sqrt (mDim);
Vector3[] dest = new Vector3[src.Length];
for (int i = 0; i < dest.Length; i++)
dest [i] = Vector3.zero;
for (int l = 0; l < bands; l++)
{
for (int mo = -l; mo <= l; mo++)
{
int outputIndex = GetIndexByLM (l, mo);
Vector3 target = Vector3.zero;
for (int mi = -l; mi <= l; mi++)
{
int inputIndex = GetIndexByLM (l,mi);
float matValue = GetValueByBand (l,mo,mi);
Vector3 source = src [inputIndex];
target += source * matValue;
}
dest [outputIndex] = target;
}
}
return dest;
}
public SHRotateMatrix(int dim)
{
mDim = dim;
mMatrix = new float[mDim][];
for (int i = 0; i < mDim; i++)
{
mMatrix [i] = new float[mDim];
for (int j = 0; j < mDim; j++)
{
mMatrix [i] [j] = 0.0f;
}
}
}
public void SetValue(int i,int j,float value)
{
mMatrix [i] [j] = value;
}
public float GetValueByBand(int l,int a,int b)
{
int centre = (l + 1) * l;
return mMatrix [centre + a] [centre + b];
}
public void SetValueByBand(int l,int a,int b,float value)
{
int centre = (l + 1) * l;
mMatrix [centre + a] [centre + b] = value;
}
private int GetIndexByLM(int l,int m)
{
return (l + 1) * l + m;
}
public int mDim;
private float[][] mMatrix;
}
这里的数学关系确实比较复杂，我也不大理解，建议直接copy即可。
4.结果
5.相关
https://github.com/TianYiJT/SphericalHarmonicLighting
http://www0.cs.ucl.ac.uk/staff/j.kautz/PRTCourse/
https://blog.****.net/baimafujinji/article/details/53869358
全局光照技术 从离线到实时渲染 第十一章
https://en.wikipedia.org/wiki/Spherical_harmonics