for (int i = 0; i<ambWidth;i++) {
for (int j = 0; j < ambHeight; j++) {
double theta = MYPI*(double)(j) / (double)ambHeight;
double phi = 2.0*MYPI*(double)(i) / (double)ambWidth;
double li = getLight(SH_Result, theta, phi);
tex_sphericaldata[(i + ambWidth*j)*ambChnnel + 0] = 255.0*li;
tex_sphericaldata[(i + ambWidth*j)*ambChnnel + 1] = 255.0*li;
tex_sphericaldata[(i + ambWidth*j)*ambChnnel + 2] = 255.0*li;
重建结果放在了tex_sphericaldata数组里，其中ambWidth是ambHeight的两倍，该结果生成一张全景图。显示全景图的程序涉及图形显示接口，太麻烦，这里就不放上了。
把图像的RGB三通道提取出来，然后求各自的球谐系数，得到下面的结果：
球谐函数的特性
SH函数有许多有趣的特性，这使得它们比我们可以选择的其他基函数更适合我们的目的。首先，SH函数不仅是正交的，而且是标准正交的，这意味着如果我们对 i 和 j 的任何一对积分
，如果i=j，计算将返回1，如果i≠j，则返回0。
SH函数也是旋转不变的，有一个定义在单位球面上的原函数：
，如果函数 g 是函数 f 的旋转副本，那么在SH投影之后：
换言之，SH投影旋转函数 g 将得到与SH投影之前将输入旋转到 f 完全相同的结果。这很令人困惑，听起来并不是什么大不了的事，但这一特性是许多其他压缩方法所不能得到的，例如JPEG的离散余弦变换编码不是平移不变的，这正是在高压缩下呈现块状外观的原因。实际上，这意味着通过使用SH函数，我们可以保证当我们设置场景动画、移动灯光或旋转模型时，照明强度不会波动、爬行（crawl）、脉冲（pulse）或有任何其他令人讨厌的瑕疵（artifacts）。
我们需要进行照明，所以一般来说，我们将对入射光进行一些描述，然后将其乘以某种表面反射率的描述（我们称之为传递函数），以得到反射光，但我们需要在入射光的整个球体上进行此操作。我们需要整合：
其中 L 是入射光，t 是传递函数。如果我们将照明函数和传递函数都投影到SH系数中，那么正交性保证函数乘积的积分与其系数的点积相同（两个函数的乘积的积分等于它们各自的球谐系数向量的点积）：
因为我们只计算了n阶的球谐函数，所以：
通过计算两个SH函数系数的点积来积分两个SH函数的乘积。
我们把球面上的积分分解成SH系数上的一个点积，只是一系列的乘法加法。这是整个过程的关键-通过将函数投影到SH空间，我们可以将球体上的积分转换为非常快速的运算。
上面的描述说明，我们想计算阴影的时候，我们有两种办法：
方法一：以前的是把照射到物体该点上的方向在 wi 的光源 Li(wi)计算其遮蔽效果，然后积分得到最终值：
方法二：现在是计算光源函数的基函数表示，然后再计算阴影传输矩阵 Mi，通过把传输矩阵 乘以 光照的SH系数，就能得到阴影光源的SH系数。
当然上面只是特性的介绍，不涉及具体求法，后面才是具体的理论。
这个点积返回一个单一的标量值，这是积分的结果，但是我们可以使用转换SH函数用在另一种技术上：
***********************************************************************************
注意，这里描述的例子，和我们目前想要做的球谐光照关系不是很大，可以忽略。
假设我们有一些我们还不知道的任意球面光源函数。我们还为曲面b(s)上的一个特定点提供了一些阴影函数，该函数描述了该点的光线是如何被遮蔽的（例如，我们上方有一个物体会阻挡来自该方向的光线），我们可以使用光线跟踪器对其进行评估。
我们需要一种方法，将入射光的系数转换成另一组被阴影函数遮住的光的系数，我们将结果称为c(s)。我们可以构造一个线性运算，使用传递矩阵将光源a(s)的SH投影直接映射到阴影光源c(s)的SH投影，而不必知道照明函数a(s)。为了建立传递矩阵M，其中矩阵的每个元素由 i 和 j 索引，计算如下：
: 计算传递矩阵元素的三重乘积。
结果是一个矩阵，我们可以使用简单的矩阵向量乘法将光源转换为阴影光源：
：将传递矩阵应用于SH系数向量
让我们尝试一个明确的例子。如果我们找到一个与SH函数完全相似的阴影函数，会发生什么？例如，如果
呢？这意味着我们将计算SH函数的三重乘积：
得到的矩阵基本上是稀疏的，25×25矩阵中非零元素的图如下所示：
图9：我们设计的三乘积转移矩阵（triple-product transfer matrix）例子中矩阵M中的非零项。
将矩阵应用于光源的SH系数，可以得到SH系数的另一个向量，其结果与在SH投影之前将光源与mask相乘的结果完全相同：
： 将示例传输矩阵应用于光源的结果。
正如你所看到的，这是一种不同的方法来记录模型上某一点的阴影，而不强迫我们做最后的积分，当我们预处理依赖于视图的光泽镜面反射曲面时，我们将使用它来获得良好的效果。
旋转球谐函数
这是整个章节中最难的一部分。
****************************************************************
首先先明白为什么需要求旋转球谐函数矩阵：
我们旋转一个三维物体的时候很简单，使用三维变换矩阵就好了。但是当我们旋转场景光源函数 f(s) 的时候，f(s) 本来可以表示成：
，如果光源转动了，那么这个时候，
而我们就需要一个系数旋转矩阵
：
就可以得到新的光源函数的表示了。（注意这里的系数
可以表示包含了阴影信息的光源函数的SH系数）
*****************************************************************
SH函数的最后一个属性是最难编程的，而且可能是大多数人都会陷入困境的地方。我们已经断言SH函数是旋转不变的，但是我们如何实际旋转SH投影函数呢？答案并不简单。首先要回答的问题是“你说的是什么样的转换方式？“你想要以欧拉角（α，β，γ）为单位的旋转吗？如果是的话，你是围绕哪个轴旋转的？”？XYZ，ZYX还是ZYZ？用四元数指定轴和角度旋转如何？把旋转推广成3×3的旋转矩阵，并有所有相关的对称冗余，怎么样？尽管Sloan, Kautz and Snyder说SH函数有“简单的旋转”，但他们并没有说明全部情况。
从正交性规则来看，SH旋转过程是一个线性运算，bands间的系数不相互作用。实际上，这意味着我们可以使用单个 n^2×n^2 旋转矩阵将SH系数向量旋转为另一个SH系数向量，并且矩阵将是块对角稀疏的，如下所示：
的确，一旦你手头上有了旋转系数，使用旋转操作就可以被视为“一个简单的计算”，但是有效地构造旋转系数远不是一件简单的事。
简而言之，让我们看看SH函数的另一种表示形式。到目前为止，我们一直用球坐标表示SH函数，但我们同样可以很容易地将其转换为（x，y，z）上的隐式函数，方法是代入球坐标到笛卡尔坐标的转换公式，并去掉项。
为此，我们提出了一组非常简单的表达式：
前几个实SH函数的笛卡尔版本:
要使用这些函数，只需在单位球体上选取一个点（x，y，z），然后通过上面的方程来计算该方向上的SH系数。用这种形式的方程来表示SH投影是可能的，但它们作为符号表达式对我们更有用。
我们可以通过构建一个矩阵来构建SH函数的旋转矩阵，其中每个元素都是通过旋转的SH样本与非旋转版本的符号积分来计算的：
这将构建一个 n^2×n^2 表达式矩阵，将SH系数的未旋转向量映射为旋转向量。例如，使用显式参数化公式：
对于前三个band，我们给出了一个绕z轴旋转的9×9矩阵：
如您所料，该矩阵扩展到更高的band，频带N使用 Nα 的正弦和余弦。
这项技术对于低阶SH函数来说很好——您只需将任何旋转分解为一系列更简单的旋转，然后重新组合结果。事实上，对于任何比二阶SH函数更大的阶数，它很快就会变成一种巨大的痛苦。
首先，我们需要多少个最小的旋转次数才能将SH函数转换成任何可能的方向？如果我们使用ZYZ公式，我们只需旋转两次（绕x轴旋转1次，再绕y轴旋转1次就可以了），其中一次我们已经有了公式！那么，如何绕y轴旋转呢？我们可以把它分解成绕x轴旋转90°，绕z轴旋转，最后绕x轴旋转-90°。很好，X轴旋转是一个固定的角度，所以我们可以将其列为一个常量浮点数组：
退一步来看这个过程的计算成本。在矩阵表示法中，我们计算：
这是4个不同的9×9矩阵乘法，加上相关的trig函数。假设矩阵矩阵乘法的代价是O（n^3），一个简单的实现将使用2916个乘法加法。我们可以利用矩阵的稀疏性，对于5阶旋转把它降低到大约340倍，但它的计算还是很耗费性能。
把旋转组合起来，再把整个运算乘法成一个大的显式表达式怎么样？对于大于1的band，这会产生一些非常可怕的trig表达式。这是前两个band的矩阵：用ZYZ-Euler角（α，β，γ）求解SH旋转矩阵前两个带的解析解
这个矩阵对调试很有用，但要使用它，我们必须将游戏引擎转换为ZYZ旋转，当旋转最终对齐时，会出现所有相关的万向节锁定问题。这不是我们都转换成四元数的原因吗？有一个诀窍可以用来阻止我们进入这个痛苦的领域，同时也可以加快计算速度。三维旋转矩阵的一个基本性质是其众多的对称性，我们可以充分利用这些对称性。给定一个普通的3×3旋转矩阵R：
我们可以使用这些恒等式直接重建ZYZ-Euler角（α，β，γ）的三角函数：
当sinβ=0时，这个公式就失效了，但在这种情况下，旋转矩阵将写为：
所以我们可以任意决定如何在α和γ之间划分总的z-旋转。强迫γ=0给出了特殊情况下的恒等式：
这两个公式都可以作为SH系数的前两个波段的权宜之计，但我们确实需要一个快速、通用的任何数量波段的解决方案。这就是我们触及当前游戏研究中未被探索的前沿。
为了解决这个问题，我们需要深入研究SH函数最初的来源。SH函数最初被设计用来在量子水平上描述单个原子中角动量的分布，这解释了为什么参数l和m是整数-它们是描述原子状态的不可分割的量子数。在使用SH函数编写程序方面经验最丰富的科学领域称为计算化学，研究人员试图在量子水平上模拟分子内部原子之间的相互作用，以更好地理解它们是如何工作的。就在1999年发表的基础研究论文描述了旋转实SH函数的新方法，这种方法比传统的Wigner D函数更有效。直到最近，这些论文的重点都集中在轴对齐的特殊情况下的显式公式（对于原子晶格的建模很好），但是我们需要一个更通用的公式来进行计算机绘图工作。
我们真的需要一组递归关系，递归函数，从带l建立带l+1的SH旋转矩阵。这样，低阶近似比高阶近似具有更低的计算复杂度。研究发现了三篇关于旋转实球谐函数递推关系的论文，它们有许多相似之处：
我已经成功地实现了Blanco的论文，但我会推荐Ivanic或Choi的论文，如果你能通过你的方法解决数学问题则它们是更有效的算法。Choi在一封个人电子邮件中说，在复数空间工作比在实数空间工作快10倍，但需要复数到实数的后处理（“复数让生活更容易！“他写道）。我正在为以后的论文实现这两种算法。有关Ivanic算法的更多信息，请参阅本文档的附录1。
（总之，太难了，所以，我们把别人写的程序抄上来就好啦~，这个工作放在最后去实现，下面先做散射光照技术）
球谐光照漫反射表面
现在我们已经介绍了SH函数的属性，并编写了一些操作它们的工具，我们终于可以使用它们为我们的模型生成一些照明。我们将通过一系列的照明技术来展示每一种照明技术是如何定义和实现的。
假设我们已经将多边形模型的描述加载到内部数据库中。我们只对使用光线跟踪器对世界进行点采样感兴趣，所以我们只需要一个顶点、法线和三角形的列表。首先，我们将顶点法线对处理成一组唯一的“照明点”，方法是复制具有多条法线的顶点，并连接在平滑曲面上共享的顶点（这可以使用STL贴图非常快速地完成）。接下来，我们将遍历模型中的所有照明点，计算每个照明点的传递函数。传递函数是一个函数，当它与一个亮度函数点乘时（即，相乘和积分）就可以表示表示该点的近似照明。
我们可以为漫反射曲面生成三种不同类型的传递函数，每种类型的传递函数的计算都会越来越复杂，因此我们将按顺序进行处理。
漫反射非阴影传输函数
回到渲染方程，我们可以将其分解为最基本的部分，即假设位于某个定向平面上的光源和曲面点，从而仅使用直接照明生成无阴影图像。
所以我们可以很好地记住，这个方程可以很好地反映光线的方向。光被均匀反射，所以光是独立于视图的，我们的视角ωo消失了。我们的BRDF只是一个简单的常量标量，可以从积分中去掉，只剩下三个元素：光源Li、简化的余弦项和线性比例因子：
（注：这里我们已经改用反照率ρx来衡量朗伯漫反射的反射率。反照率是发射的辐射强度与辐照度的比值，对于漫反射表面，它会收缩为πρx，因此我们只需要指定一个标量ρx，它的值在[0,1]范围内。就像RGB颜色一样。）
将光源从传递函数中分离出来（对于漫反射非阴影传输，我们称之为MDU），我们得到了用SH函数近似的函数。它只是一个几何项，法线和光源之间夹角的余弦，固定在零：
为了计算这个函数，我们将把预先计算好的光线和SH系数列出来，并遍历它们，在每个采样点上用曲面法线打点，看看它是否在上半球内部。
除了使用任意复杂的SH区域光源外，这将产生与普通点积照明非常相似的图片。使用5阶漫反射非阴影SH传递函数（25个系数）进行渲染。实际上，光照仅使用几何项：
自己写程序的话，首先生成几何项的球谐系数。这里为了方便起见，不给每个面生成，而是直接给屏幕中采样到的结果来生成：
注意：我写的程序里，生成的图像的长宽都是 ThreadNum
这里的渲染框架用的是光线追踪三部曲中的程序：《Ray Tracing in One Weekend 》等三本小教程
不过如果要构建球谐光照的程序的话，其实搞完第一本书就好了。后面的包围盒、加速结构之类的暂时其实不需要。可以根据自己的兴趣，在里面实现三角面片的渲染。
void GeometryInit() {
SH_Geometry = new double*[ThreadNum*ThreadNum];
for (int i = 0; i < ThreadNum; i++) {
for (int j = 0; j < ThreadNum; j++) {
SH_Geometry[i + ThreadNum*j] = new double[3 * n_bands*n_bands];
float u = (float(i) + 0.5) / float(ThreadNum);//random()
float v = (float(j) + 0.5) / float(ThreadNum);//random()
Ray r = cam.get_ray(u, v);
//i-width j-height
GeometryGenerate(r, world, SH_Geometry[i + ThreadNum*j]);
对于采样到的屏幕上的每个点，都计算其球谐系数：
double ** SH_Geometry;
int red_offset = 0;
int green_offset = 0 + n_bands*n_bands;
int blue_offset = 0 + 2*n_bands*n_bands;
void GeometrySpherical(double *GeometryPara, const SHSample samples[], Vec3d normal,Vec3f albedo) {
int n_samples = sqrt_n_samples*sqrt_n_samples;
int n_coeff = n_bands*n_bands;
// for each sample
for (int i = 0; i<n_samples; ++i) {
// calculate cosine term for this sample
double H = dot(samples[i].vec, normal);
if (H > 0.0) {
// ray inside upper hemisphere so...
// SH project over all bands into the sum vector
for (int j = 0; j<n_coeff; ++j) {
double value = H * samples[i].coeff[j];
GeometryPara[j + red_offset] += albedo.x * value;
GeometryPara[j + green_offset] += albedo.y * value;
GeometryPara[j + blue_offset] += albedo.z * value;
else {
// ray not in upper hemisphere
// divide the result by probability / number of samples
double factor = 1.0 / (double)n_samples;//注意都是整形，一定要变成浮点数。
for (int i = 0; i<3 * n_coeff; ++i) {
GeometryPara[i] = GeometryPara[i] * factor;
void GeometryGenerate(const Ray&r, hitable *world, double *GeometryPara) {
int n_coeff = n_bands*n_bands;
for (int i = 0; i<3 * n_coeff; ++i) {
GeometryPara[i] = 0;
hit_record hrec;
if (world->hit(r, 0.0001, MAXFLOAT, hrec)) {
scatter_record srec;
Vec3f emitted = hrec.mat_ptr->emitted(r, hrec, hrec.texU, hrec.texV, hrec.p);
hrec.mat_ptr->scatter(r, hrec, srec);
Vec3d Normal(hrec.normal.x, hrec.normal.y, hrec.normal.z);
GeometrySpherical(GeometryPara, SHSampleList, Normal, srec.attenuation);
else {
最终渲染显示界面的话就简单了：
for (int i = 0; i < ThreadNum; i++) {
for (int j = 0; j < ThreadNum; j++) {
Vec3d col(0.0,0.0,0.0);
int offset = (i + ThreadNum * j);
for (int k = 0; k<n_coeff; ++k) {
col.x += SH_Geometry[offset][k + red_offset] * SH_ResultR[k];
col.y += SH_Geometry[offset][k + green_offset] * SH_ResultG[k];
col.z += SH_Geometry[offset][k + blue_offset] * SH_ResultB[k];
clamp(col.x, 0.0, 1.0);
clamp(col.y, 0.0, 1.0);
clamp(col.z, 0.0, 1.0);
myHostFrameBuffer.finalBuffer_f[offset * A_PIXEL_BYTES_NUM + 0] += col.x;
myHostFrameBuffer.finalBuffer_f[offset * A_PIXEL_BYTES_NUM + 1] += col.y;
myHostFrameBuffer.finalBuffer_f[offset * A_PIXEL_BYTES_NUM + 2] += col.z;
最终渲染结果如下：
在屏幕空间做球谐系数的坏处是只求一次系数的话不抗锯齿。所以我求了6次，得到平均的球谐系数：
void GeometryInit() {
int n_coeff = n_bands*n_bands;
SH_Geometry = new double*[ThreadNum*ThreadNum];
for (int i = 0; i < ThreadNum; i++) {
for (int j = 0; j < ThreadNum; j++) {
SH_Geometry[i + ThreadNum*j] = new double[3 * n_bands*n_bands];
for (int k = 0; k<3 * n_coeff; ++k) {
SH_Geometry[i + ThreadNum*j][k] = 0;
for (int count = 0; count < 6; count++) {
for (int i = 0; i < ThreadNum; i++) {
for (int j = 0; j < ThreadNum; j++) {
double* tempSHGeo = new double[3 * n_bands*n_bands];
float u = (float(i) + random()) / float(ThreadNum);//
float v = (float(j) + random()) / float(ThreadNum);//
Ray r = cam.get_ray(u, v);
//i-width j-height
GeometryGenerate(r, world, tempSHGeo);
for (int k = 0; k<3 * n_coeff; ++k) {
SH_Geometry[i + ThreadNum*j][k] += tempSHGeo[k]/(double)6;
我换用一下书里给的光照，就是这个：
得到的结果是：
可以看到，在这个视角下，光源在左边的方向。
我采用了9阶band，加15次渲染几何球谐系数取平均（大概用了十几分钟来渲染），得到结果：
漫反射阴影传输函数
在简化的渲染方程中加入可见项，我们可以将自阴影添加到照明模型中，我们的照明开始变得非常有趣。
只有这个微小的变化，物体表面上的点不再被假设为位于一个无限平面上，而是与周围的几何体相互作用，并且可以访问其入射区域光源的视野。漫反射阴影照明 M DS 的结果传递函数为：
这一效果，有时被称为遮挡环境，是传统CG和全局照明图像之间最强大的区别。计算传递函数需要我们从当前点通过多边形数据库跟踪光线，以找到任何命中-请注意，我们不需要来自命中的任何几何信息，只需要它发生的布尔值。
// for each sample
for (int i = 0; i<n_samples; ++i) {
// calculate cosine term for this sample
Hs = DotProduct(sample[i].vec, normal);
if (Hs > 0.0) {
// ray inside upper hemisphere...
if (!self_shadow(pos, sample[i].vec)) {
// ray hits nothing, add in its the contribution:
for (int j = 0; j<n_coeff; ++j) {
value = Hs * sample[i].coeff[j];
result[j + red_offset] += albedo_red * value;
result[j + green_offset] += albedo_green * value;
result[j + blue_offset] += albedo_blue * value;
else {
// ray hits self...
else {
// ray not in upper hemisphere
// divide the result by number of samples
double factor = area / n_samples;
for (i = 0; i<3 * n_coeff; ++i) {
coeff[i] = result[i] * factor;
self_shadow()函数取决于光线跟踪器，但是对于SH投影过程，有一些事情需要注意，这使得代码更易于编写。
首先，如果决定使用像体素栅格、层次边界盒树或BSP树之类的加速数据结构，请记住每个光线原点都将位于对象的边界框内，因此不需要初始光线盒交点来找到光线跟踪的起点。只需将起始点散列到数据结构中，然后从那里开始。第二，self_shadow的重点是找到遮挡，所以我们希望我们的阴影触角击中物体，我们将从模型上的一个顶点向各个方向发射光线。
问题是，与顶点相邻的多边形总是与许多光线共面。对包含光线原点作为顶点的多边形进行测试的任何光线最多将返回原点处的命中（需要进行epsilon测试以排除距离原点太近的命中），或者最坏的情况下返回沿光线的某个位置的命中。如果不小心，这可能会导致对象上出现不正确的阴影，因此，如果可能，请在加载模型时保留一些面邻接信息，并从阴影测试中排除共享当前顶点的多边形。
当从多边形模型上的顶点跟踪光线时，光线多边形问题的示例。
第三个问题是单面多边形测试。为了使光线多边形测试正确返回遮挡，我们不能使用单面多边形测试。排除光线多边形测试，因为如果光线是从模型内开始的（上图中的A和C相交），则多边形与光线“背离”的情况将不起作用，这与许多阴影光线的情况相同。这一限制还意味着，一般来说，我们必须有多个对象——具有完全封闭的皮肤的对象。如果相交处没有顶点，具有自相交曲面的对象将具有不正确的阴影（在使用Gouraud着色重建后）。当照亮单一厚度的物体时要小心，确保任何阴影触角都会被其他物体正确遮挡（因为下面的插图在后壁上没有正确地进行！）
阴影光线还可以在找到遮挡多边形后立即返回交点（交点B），因为相交的顺序并不重要，只需要存在一个交点。
使用5阶漫反射阴影SH传递函数的渲染。注意来自恒定半球光源的柔和阴影 :
为了简单我就不写上面的优化方法了，直接搞个简单粗暴的版本：
void GeometrySpherical(double *GeometryPara, const SHSample samples[], Vec3d normal,Vec3f& albedo, Vec3f& scatteredPoint, hitable *world) {
int n_samples = sqrt_n_samples*sqrt_n_samples;
int n_coeff = n_bands*n_bands;
// for each sample
for (int i = 0; i<n_samples; ++i) {
// calculate cosine term for this sample
double H = dot(samples[i].vec, normal);
if (H > 0.0) {
// ray inside upper hemisphere so...
// SH project over all bands into the sum vector
Vec3f scatteredVec(samples[i].vec.x, samples[i].vec.y, samples[i].vec.z);
Ray scattered = Ray(scatteredPoint, scatteredVec, 0.0);
hit_record hrec;
if (!world->hit(scattered, 0.0001, MAXFLOAT, hrec)) {
for (int j = 0; j<n_coeff; ++j) {
double value = H * samples[i].coeff[j];
GeometryPara[j + red_offset] += albedo.x * value;
GeometryPara[j + green_offset] += albedo.y * value;
GeometryPara[j + blue_offset] += albedo.z * value;
else {
// ray not in upper hemisphere
// divide the result by probability / number of samples
double factor = 1.0 / (double)n_samples;//注意都是整形，一定要变成浮点数。
for (int i = 0; i<3 * n_coeff; ++i) {
GeometryPara[i] = GeometryPara[i] * factor;
void GeometryGenerate(const Ray&r, hitable *world, double *GeometryPara) {
int n_coeff = n_bands*n_bands;
for (int i = 0; i<3 * n_coeff; ++i) {
GeometryPara[i] = 0;
hit_record hrec;
if (world->hit(r, 0.0001, MAXFLOAT, hrec)) {
scatter_record srec;
hrec.mat_ptr->scatter(r, hrec, srec);
Vec3d Normal(hrec.normal.x, hrec.normal.y, hrec.normal.z);
GeometrySpherical(GeometryPara, SHSampleList, Normal, srec.attenuation, hrec.p, world);
给三角模型求交计算太慢了，我运行了半个小时都没出结果，醉了，直接几个最简单的模型来渲染吧：
因为这个正方体就是凌空放的，所以阴影和立方体之间有一点缺口。（和光栅引擎的阴影映射加bias纠正阴影深度的效果很像，注意区分一下。。。。）
因为没有优化，所以就这几个破模型我渲染了十几分钟。。。。
我把模型改了改,不让它浮空了，然后又换了个场景光源：
漫反射相互反射传输函数
最后一种漫反射照明方法是最引人注目的。渲染方程中有趣的一点是，它递归地添加不是直接来自光源的光，而是作为从模型上某个点可见的其他多边形的二次反射光。我们可以用积分来表示：
相互反射光传递函数很难用数学方法简洁地描述，也没有真正的照明效果，但是生成它的算法更容易描述。有四个步骤：
a） 对于模型上的每个着色点x，计算该点的直接照明传递函数（即上一节中的漫反射阴影照明）。
b） 接下来，从当前点发射一条光线，直到其中一条光线击中对象上的另一个三角形。使用命中的重心坐标在三角形的每个角上线性插值SH函数。此传递函数是指反射回着色点的光线量。
c） 将反射光乘以光线与曲面法线在x处的点积（矢量标量乘，根据余弦项缩放此相互反射光的反射量），并将其相加为一个空的SH向量。所有光线投射完毕后，像往常一样，将累积值除以采样数和蒙特卡罗加权项。
d） 一旦所有的着色点都被计算出来，这个新的SH向量集就是一个漫反射光的反弹，而且只有相互反射光。对于其他反弹，请使用此新值集作为每个三角形顶点处的起始灯光强度重复光线跟踪。重复这个动作，直到没有能量传输，或者直到你达到N次反弹。最后，将所有反弹加上直接照明相加为一个SH向量列表。
一个漫反射自转移的样本，向下的平面A发射一条光线。
我们可以看到，渲染出来的地面有点发红，这是多次反射的结果。
鉴于这个程序写出来我估计得跑一整天才能出结果，人生这么短暂，我又不是球谐渲染专家，也不从事球谐渲染研究，而且最近心情不好，不想等这么久，总之我选择跳过。
渲染SH漫反射表面
现在我们已经为每个顶点设置了一组SH系数，我们如何使用当前图形硬件来构建渲染器，以提供实时帧速率？回到SH函数的性质，SH照明的基本计算是SH投影光源与SH传递函数之间的点积：
这个计算将为我们提供一个顶点的单通道光强度，并且我们通过使用Gouraud着色填充顶点之间的间隙来近似对象的完整解。请记住，所有的SH计算都发生在对象空间中，因此如果要在世界空间中重定向对象或旋转照明功能，则需要先将灯光旋转到对象空间中。
在这一点上，我们可以做出一些选择。如果我们假设照明是白色光源（红色、绿色和蓝色相同），则每个顶点的计算值为：
for(int j=0; j<n_coeff; ++j) {
vertex[i].red += light[j] * vertex[i].sh_red[j];
vertex[i].green += light[j] * vertex[i].sh_green[j];
vertex[i].blue += light[j] * vertex[i].sh_blue[j];
对于每个颜色通道，彩色光源将具有不同的SH函数，进行计算：
for(int j=0; j<n_coeff; ++j) {
vertex[i].red += light_red[j] * vertex[i].sh_red[j];
vertex[i].green += light_green[j]* vertex[i].sh_green[j];
vertex[i].blue += light_blue[j] * vertex[i].sh_blue[j];
请注意，在这两种计算中都没有提到曲面法线-它已经隐式编码在SH系数中，就像环境项一样，因此我们使用内置的阻挡环境光来获得漫反射着色。
一个优化是要注意到，如果模型不太需要强调硬阴影（比如人脸），余弦项可以很好地通过前两个阶的SH函数来近似，从而为我们提供了一个4系数的下降来代替正常的漫反射着色项，它为我们提供了软的自阴影，而不需要为每个顶点增加额外的乘法没有额外的存储空间。这个4系数系统也可以很好地转化为4路SIMD操作，使用与3个平行光源加环境矩阵完全相同的代码。
另一个选择是使用单通道白光源，只编码自阴影，忽略自转移。最终的照明计算只需要一个传递函数和每个顶点一种颜色，并且可以动态地重新着色：
for(int j=0; j<n_coeff; ++j) {
vertex[i].red += k_red * light[j] * vertex[i].sh[j];
vertex[i].green += k_green * light[j] * vertex[i].sh[j];
vertex[i].blue += k_blue * light[j] * vertex[i].sh[j];
我们可以混合和匹配SH照明和普通照明，因为灯光是线性相加的，所以我们可以将SH照明用于着色器的环境光和漫反射部分，并在顶部添加一个镜面反射项。这正是我们在本文档开始时对汽车图像所做的，首先用SH光和alpha混合在顶部的预过滤镜面反射环境贴图来渲染车身。
到目前为止，我们已经在一个物体上使用了相同的无限光源，但是有一种方法可以模拟照亮模型的局部光源。这只在我们只有自阴影和没有自转移的情况下才起作用，因为它打破了模型中低照明方差的假设。为了进行局部照明，我们从模型表面上几个间隔良好的预计算点构造局部照明环境的少量样本（例如6或8），计算光照函数，就像模型不存在一样。如果每个球形样品附近都有发射器，则会在稍微不同的位置看到局部光源。我们通过计算每个顶点的照明采样的加权和，以顶点到灯光采样点的距离为权重，重建正确的局部照明。斯隆、考茨和斯奈德的SH论文中的例子使用了加权方案：
在顶点 i 和照明点 j 之间，其中n=10，在Max Planck的示例图像中，强调局部照明的效果。其他研究人员发现了一种更快、更宽松的方法。取N个照明采样位置，构造一个非常低的多边形三角形多面体。下一步，根据最接近的多边形将模型分成多个顶点组。为了预先计算光在一组顶点上的分布，我们将顶点投影到照明三角形上，并计算该投影点的二维重心坐标。然后我们用三个光源的权重。这种技术的另一个优点是可以在硬件上进行计算。
我们只触及了重建SH照明技术的表面，所以试试看你能不能想出一些更好更灵活的效果。目前的圣杯是将SH照明从静态模型技术扩展到可以在动画模型上工作的技术。混合SH向量的方法与我们对关节使用混合矩阵的方式对余弦项效果很好，但是阴影函数V太随机和不可预测，不能表示为几个快照的线性组合。这个问题还需要做更多的工作。
1955年，负责定义照明和色彩科学的国际标准机构CIE发表了一篇论文，定义了判断建筑设计的三种标准参考光源。该标准可从CIE获得，价格约为10欧元，并且是版权，所以我只能给你从网上窃取的单频道版本。从CIE阴天模型开始，该模型提供均匀阴天的照明功能：
CIE晴空模型有点复杂，因为我们必须考虑到太阳的位置以及它接近地平线时的不同散射效应。
相当令人印象深刻的等式。产生的光源非常逼真。
CIE部分云层模式是介于两者之间的中间地带，比晴空模式简单一些。
晴空模式的定义仍然有效。
我们可以使用的另一种照明模型是保罗·德贝维奇（Paul Debevec）制造的高动态范围光探头。HDR图像是使用相机校准和多次曝光的特殊过程从真实世界捕获的浮点RGB值的编码位图，该过程对场景中的全部光能进行编码。使用免费提供的程序HDRShop，我们可以生成HDR场景的角度映射投影，并将它们保存为float RGB像素数组。
因为是生成光源嘛，所以其实没啥意思，我们一般都用全景图来渲染，比如
http://www.hdrlabs.com/sibl/archive.html 这个网站有一堆全景图，我选了个这个图：
渲染结果为：
可恶，照明几乎全都是用的中间的黄色天际线（因为图中正方形的法线方向y分量是0），导致整个界面看起来黄黄的。。。。
为了试个高级效果，我决定加上bunny这个小兔子，把渲染图改改：
好了，等上一晚上看看明天能不能出结果吧~
终于出结果了：
高级的SH技术
光泽表面的SH照明将包含在本文的扩展版本中，很快可从以下网站获得：
http://research.scea.com/
希望结合注释和示例代码：
•解开坐标空间的混合。
•使用图形硬件计算SH传输函数。
•扩展混合材料模型的预处理器。
•扩展静态发射表面的预处理器。
•向预处理器添加反射焦散。
•漫反射体照明，可实时重新照亮云层。
•投影传输，用于采样对象周围的阴影。
•实时半透明。
希望我们已经展示了SH照明如何使用一点额外的处理能力，在静态和动态模型上生成非常逼真的图像。我们已经展示了如何使用交互反射对静态模型进行预处理，如何生成漫反射阴影的传递函数，以及它如何影响模型的设计方式。SH照明可以作为静态对象的普通渲染器的漫反射和环境光的替代品，我们已经介绍了实时更新全局和局部照明的一系列选项。
我希望你能为你的下一个项目找到一个SH照明的用途，这样我们就可以把游戏图形带到现实主义的下一个阶段。我很期待看到人们如何使用SH函数来完成需要编码和操作球形函数的意外任务。记得写下来分享你的想法！
参考文献：
《Green R. Spherical harmonic lighting: The gritty details[C]//Archives of the Game Developers Conference. 2003, 56: 4.》
代码附录一
实现球谐投影和重建，在main中打印出来球谐系数。（虽然程序大部分不是我自己写的，但是把它们组合在一起，能够直接测试成功，对接下来的学习会有很大的帮助！）
需要自己实现一个Vec3d类，表示三维double型向量。以及生成0-1随机数的程序random()，以及 MYPI: 3.14159265359
除此之外不需要包含什么头文件。
//注意0的阶乘是1
// Number of precomputed factorials and double-factorials that can be
// returned in constant time.
const int kCacheSize = 16;
double Factorial(int x) {
static const double factorial_cache[kCacheSize] = { 1, 1, 2, 6, 24, 120, 720, 5040,
40320, 362880, 3628800, 39916800,479001600, 6227020800,87178291200, 1307674368000 };
if (x < kCacheSize) {
return factorial_cache[x];
else {
double s = factorial_cache[kCacheSize - 1];
for (int n = kCacheSize; n <= x; n++) s *= n;
return s;
double DoubleFactorial(int x) {
static const double dbl_factorial_cache[kCacheSize] = { 1, 1, 2, 3, 8, 15, 48, 105,
384, 945, 3840, 10395, 46080,
135135, 645120, 2027025 };
if (x < kCacheSize) {
return dbl_factorial_cache[x];
else {
double s = dbl_factorial_cache[kCacheSize - (x % 2 == 0 ? 2 : 1)];
double n = x;
while (n >= kCacheSize) {
s *= n;
n -= 2.0;
return s;
const int n_bands = 16;
struct SHSample {
Vec3d sph;
Vec3d vec;
double *coeff;
const int sqrt_n_samples = 100;
SHSample SHSampleList[sqrt_n_samples*sqrt_n_samples];
double P(int l, int m, double x)
// evaluate an Associated Legendre Polynomial P(l,m,x) at x
double pmm = 1.0;
if (m>0) {
double somx2 = sqrt((1.0 - x)*(1.0 + x));
double fact = 1.0;
for (int i = 1; i <= m; i++) {
pmm *= (-fact) * somx2;
fact += 2.0;
if (l == m) return pmm;
double pmmp1 = x * (2.0*m + 1.0) * pmm;
if (l == m + 1) return pmmp1;
double pll = 0.0;
for (int ll = m + 2; ll <= l; ++ll) {
pll = ((2.0*ll - 1.0)*x*pmmp1 - (ll + m - 1.0)*pmm) / (ll - m);
pmm = pmmp1;
pmmp1 = pll;
return pll;
double K(int l, int m)
// renormalisation constant for SH function
double temp = ((2.0*l + 1.0)*Factorial(l - m)) / (4.0*MYPI*Factorial(l + m));
return sqrt(temp);
double SH(int l, int m, double theta, double phi)
// return a point sample of a Spherical Harmonic basis function
// l is the band, range [0..N]
// m in the range [-l..l]
// theta in the range [0..Pi]
// phi in the range [0..2*Pi]
const double sqrt2 = sqrt(2.0);
if (m == 0) return K(l, 0)*P(l, m, cos(theta));
else if (m>0) return sqrt2*K(l, m)*cos(m*phi)*P(l, m, cos(theta));
else return sqrt2*K(l, -m)*sin(-m*phi)*P(l, -m, cos(theta));
void clamp(double& data, double min, double max) {
if (data < min)data = min;
if (data > max)data = max;
double getLight(double result[], double theta, double phi) {
double data = 0.0;
for (int l = 0; l<n_bands; ++l) {
for (int m = -l; m <= l; ++m) {
int index = l*(l + 1) + m;
//注意m是从-1开始的
data += result[index]*SH(l, m, theta, phi);
clamp(data,0.0,1.0);
return data;
void SH_setup_spherical_samples(SHSample samples[], int sqrt_n_samples)
// fill an N*N*2 array with uniformly distributed
// samples across the sphere using jittered stratification
int i = 0; // array index
double oneoverN = 1.0 / sqrt_n_samples;
for (int a = 0; a<sqrt_n_samples; a++) {
for (int b = 0; b<sqrt_n_samples; b++) {
// generate unbiased distribution of spherical coords
double x = (a + random()) * oneoverN; // do not reuse results
double y = (b + random()) * oneoverN; // each sample must be random
double theta = 2.0 * acos(sqrt(1.0 - x));
double phi = 2.0 * MYPI * y;
samples[i].sph = Vec3d(theta, phi, 1.0);
// convert spherical coords to unit vector
Vec3d vec(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta));
samples[i].vec = vec;
samples[i].coeff = new double[n_bands*n_bands];
// precompute all SH coefficients for this sample
for (int l = 0; l<n_bands; ++l) {
for (int m = -l; m <= l; ++m) {
int index = l*(l + 1) + m;
//注意m是从-1开始的
samples[i].coeff[index] = SH(l, m, theta, phi);
double SH_Result[n_bands*n_bands];
//theta 0-pi phi 0-2pi
typedef double(*SH_polar_fn)(double theta, double phi);
void SH_project_polar_function(SH_polar_fn fn, const SHSample samples[], double result[])
int n_coeff = n_bands*n_bands;
const double weight = 4.0*MYPI;
int n_samples = sqrt_n_samples*sqrt_n_samples;
// for each sample
for (int i = 0; i<n_samples; ++i) {
double theta = samples[i].sph.x;
double phi = samples[i].sph.y;
for (int n = 0; n<n_coeff; ++n) {
result[n] += fn(theta, phi) * samples[i].coeff[n];
// divide the result by weight and number of samples
double factor = weight / n_samples;
for (int i = 0; i<n_coeff; ++i) {
result[i] = result[i] * factor;
/****************** 测试区 *****************************/
double getLightIntensity(double theta, double phi) {
double intensity = std::max(0.0, 5 * cos(theta) - 4) + std::max(0.0, -4 * sin(theta - MYPI)*cos(phi - 2.5) - 3);
return intensity;
int main(void){
SH_setup_spherical_samples(SHSampleList, sqrt_n_samples);
SH_project_polar_function(getLightIntensity, SHSampleList, SH_Result);
/*把 SH_Result 的结果打印出来，如果是和
0.398787
-0.210511 0.286563 0.281813
-0.315004 -1.0256e-5 0.131397 9.08118e-5 0.0931611
-0.249618 -8.40392e-5 0.123395 0.303834 -0.165057 3.75637e-5 -0.092273
差不多，说明你运行的结果就是对的。
游戏中的光照可以被简单分为两个部分，直接光和间接光。
直接光中我们研究各种不同的BRDF材质，甚至BSDF，BSSSRDF等等。这些模型据有很不错的表现力，足够我们区分金属，皮肤，木头等等不同物体的着色表现。
但这并不能满足我们，因为光并不是那么简单，光会被反射，会被折射，会被透射，会被吸收，所以物体的受光情况同时又由这个场景的其他物体决定，这部分光照同时拥有着更加富强的表现力，...
一、探针基于球谐函数的全局光照
球谐光照是基于预计算辐射度传输理论实现的一种实时渲染技术。预计算辐射度传输技术能够重现在区域面光源照射下的全局照明效果。这种技术通过在运行前对场景中光线的相互作用进行预计算，计算每个场景中每个物体表面点的光照信息，然后用球谐函数对这些预计算的光照信息数据进行编码，在运行时读取数据进行解码，重现光照效果。
球谐光照使用新的光照方程来代替传统的光照方程，并将这些新方...
2. 球谐变换
类比傅里叶变换[采用定义在圆上的三角函数]，球谐变换采用定义在球面上的一组球谐基函数。[Spherical Harmonic Lighting: The Gritty Details]有这些函数的介绍。
现在只需要知道
球谐光照（Spherical Harmonics Lighting）
文章目录
球谐光照（Spherical Harmonics Lighting）一、前言二、
球谐函数2.1 基函数2.2 投影与重建2.3 应用三、漫反射环境光3.1 IrradianceMap3.2 数学推导3.3 实践参考博文
在学习图形渲染的过程中，一直對
球谐函数（
球谐光照）有一点了解，但没有亲手实现过。
关于
球谐函数的数学概念是比较复杂的，笔者也实在无法完全理解，这篇文章只能从做东西的角度来说，公式是如何和代码进行对应的。
一、球谐光照
球谐光照(Spherical Harmonic Lighting)就是基于球面调和(SH, Spherical Harmonics)这个数学工具的一种光照/着色算法。球谐光照也是一种用让我们捕捉(capture)光照、随后进行重新光照(relight)、实时展现全局光照(Global Illumination)风格的区域光源(area light sources)与软阴影的一种技术。
如果你可以理解一维函数的傅立叶变换(Fourier Transform)，那么其实你也很容易理解在图形学领域里
文章目录使用球谐函数的本质蒙特卡洛积分正交基函数球谐(Spherical Harmonics)球谐投影(SH Projection)SH函数的特性(Properties of SH Functions)Rotating Spherical HarmonicsSH光照表面漫反射(SH Lighting Diffuse Surfaces)Envmap的例子说明参考文献
使用球谐函数的本质
使用球谐函数的本质就是把一个点上的空间光照信息投影到球谐基向量上，然后再在需要使用的时候，使用投影得到的球谐系数来还原逼近
OpenGL 半球光照原理是一种基于环境光照的光照模型，它使用半球体来模拟环境光照射到物体表面的效果。半球体通常被分成两个部分，一个上半球和一个下半球，上半球表示天空光照，下半球表示地面光照。在这种光照模型中，每个顶点都会计算其对应的法向量，然后使用法向量去采样半球体上的纹理，从而得到该点的环境光颜色。
下面是一个简单的OpenGL半球光照的Shader代码：
顶点着色器：
#version 330 core
layout (location = 0) in vec3 aPos;
layout (location = 1) in vec3 aNormal;
out vec3 vPos;
out vec3 vNormal;
uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;
void main()
gl_Position = projection * view * model * vec4(aPos, 1.0);
vPos = vec3(model * vec4(aPos, 1.0));
vNormal = mat3(transpose(inverse(model))) * aNormal;
片段着色器：
#version 330 core
out vec4 FragColor;
in vec3 vPos;
in vec3 vNormal;
uniform vec3 lightColor;
uniform vec3 objectColor;
uniform vec3 lightPos;
uniform vec3 viewPos;
uniform samplerCube skybox;
void main()
vec3 ambient = 0.2 * lightColor;
vec3 norm = normalize(vNormal);
vec3 lightDir = normalize(lightPos - vPos);
float diff = max(dot(norm, lightDir), 0.0);
vec3 diffuse = diff * lightColor;
vec3 viewDir = normalize(viewPos - vPos);
vec3 reflectDir = reflect(-lightDir, norm);
float spec = pow(max(dot(viewDir, reflectDir), 0.0), 32.0);
vec3 specular = spec * lightColor;
vec3 envColor = texture(skybox, normalize(vPos)).rgb;
vec3 color = vec3(0.0);
color += ambient * envColor;
color += diffuse * objectColor * envColor;
color += specular * envColor;
FragColor = vec4(color, 1.0);
在这个Shader中，我们使用了天空盒纹理来模拟环境光照。通过采样天空盒纹理，我们可以得到当前顶点的环境光颜色。然后我们计算漫反射光和镜面反射光的强度，最终得到该点的最终颜色。
Input type (torch.cuda.DoubleTensor) and weight type (torch.cuda.FloatTensor) should be the same
贝_探索版:
PBRT——零基础到完全吃透系列
Dezeming:
VSCode与Latex环境的搭建（最简洁，最省事，最舒服的方案，不用搞一堆乱七八糟的配置）
quaer:
错误 MSB8066 C:\Program Files (x86)\Microsoft Visual Studio\2019\Professional\MSBuild\Microsoft\VC\v16
python激活python38
什么是CCS Concepts
就是用的复指数基来进行积分的。
作为迭代循环的初始状态，它把我们要求的目标多项式从0阶提升到m阶。
，如果
就正好结束。如果 m<l ，就不断把 m 提升来达到 l：
左图是光源的颜色，可知上面和靠右边的部分很亮。
，因此权重就是 