【Siggraph 2001】An Efficient Repr

作者: 离原春草 | 来源:发表于2021-08-13 10:31 被阅读0次

这里给出的是Siggraph 2001年给出的一项使用球谐来表示Irradiance Environment Map的技术，虽然年代已经十分久远了，但是这项技术至今仍在众多的引擎中发挥着重要作用，本文的相关论文以及PPT资源在参考文献中可以找到。

个人总结

这篇文章提出了一种使用SH来表达天光Irradiance Map的算法，通过这个算法，Irradiance Map可以表示成9个系数（RGB分量各9个）的SH矩阵，而光照的计算则从此前的采样变成了法线与矩阵的乘法，不但大大降低了内存消耗，同时还减少了贴图采样导致的时间消耗，算是一个十分经典的简化做法了。

正文

如上图所示，环境光对光照效果有着十分重要的意义，当然，如果不做简化处理，渲染的时候添加环境光，其消耗也是十分之高昂的。

先来看下前人是如何实现环境光反射的：

Blinn & Newell在1976年的时候提出使用reflection map的方式来实现环境光反射，从上面的截图判断，reflection map应该是一张2D贴图（这张贴图覆盖了各个方位的远景光照输入，跟cubemap也比较相似了）。

与之类似的技术有environment map，这是一张烘焙了各个方向光照强度的cubemap，根据任意点的法线与视线，可以很容易的找到反射点的光照数据，因此可以很方便的实现镜面反射效果。

Environment Map

Irradiance Map

对于漫反射而言，我们需要采集的就不是一个反射方向的光照输入，而是法线方向上半球面的所有输入，而这个输入由于只跟法线有关，因此可以放在离线完成，计算的结果可以用Irradiance Map表示，在运行时只需要根据法线对这个cubemap进行采样即可得到对应的diffuse输入数据。

如上图所示，简单来说，漫反射可以看成是反射响应函数（基色贴图）与输入照度（Irradiance）之间的混合，基色贴图是可以直接获取的，但是输入光照度的计算就比较复杂了。

这里的一个问题是，为什么diffuse反射计算用到的是Irradiance而非Radiance，实际上，前者是单位面积下的光辐射通量，后者是单位面积单位固体角下的光辐射通量，也就是两者的区别在于前者是没有方向性的，而后者是有方向性的，因此这里推测，在Lambert Surface上面，在同样的输入下，各个方向上的输出光强是一致的，也就是说是没有方向性的，因此这里会使用输入Irradiance Map而非Radiance Map（我们计算漫反射，需要考虑的是各个不同方向的输入光照作用下的总的反射结果，而非仅仅单一方向的输入光的反射结果）。

此前照度计算有如上所示的众多方案，而照度实际上在物件表面的变化是十分缓慢的，因此我们在实际使用中可以考虑使用低分辨率贴图来承载照度数据。

除了前面的测量解之外，照度计算还可以通过解析的方式完成，比如这里说到Radiance/Irradiance数据可以分解成多个basis function（球谐函数）的，各个basis function按照一定的系数进行累加就能得到原始数据的近似与模拟。

对于某个位置的Irradiance数据，可以只需要9个系数就能完成比较精确的表达，而这些系数可以通过二次多项式计算得到。

这个方法相对于以往的方案有几点创新：

提出了频域的分析方法
给出了Irradiance的高效计算方案
可以用于实现程序化的渲染，不需要光照贴图的参与
给出了一种用于表达光照的全新方式

传统的Irradiance计算方式是通过为每个像素进行一次半球积分来计算得到，其结果相当于将高频数据过滤掉了，只保留了低频信息，而按照这种思路，我们完全可以将这个工作放到频域完成，其效率更高。

球谐函数是一组相互正交的函数，这些函数叫做基函数，也就是前面说的basis function，如上图所示，这里只考虑前三阶，第0阶是1，第一阶对应三个坐标轴的一阶函数，第二阶则对应三个坐标轴的二阶函数。

基函数两两相乘之后的积分结果为0，而与本身相乘之后的积分结果为1，这是球谐函数的基本性质。

光照L与照度E可以通过与球谐函数相乘并积分来得到这些信号在对应基函数上的投影系数，从而通过这些投影系数与对应基函数相乘累加就能实现对应信号的近似模拟。

这里给了一个在Lambertian Surface上从光辐射率（L）的球谐系数到光照度（E）球谐系数的转换公式（不知道怎么推导出来的，理论上 $L=\frac{dE}{d\omega}$ ，只有当反射点上各个角度接受到的光辐射率都相同，这个公式才会成立，由于这里说的是Lambertian Surface，也就意味着光照打在这个物件上之后，朝着各个方向的散射光强是相同的，而由于我们在计算物件光照的时候，需要的L跟E都是出射，因此这个公式的条件就满足了），这个转换是通过乘上一个换算函数 $A_l$ 实现的，而这个换算函数实际上就是一个低通滤波器，保留低频部分，移除高频部分，也就是说，这里通过一个低通滤波乘法就完成了前面的半球积分。

这里给出了不同阶数的SH模拟的误差数据，可以看到阶数越高，精度越高，而到了3阶，其误差已经降低到1%，基本上已经肉眼很难分辨了，这也是为什么我们在使用SH进行光照表达时，不会超过3阶的原因。

另外这里我们看到，在本文中的基函数（多项式）是不带任何前置系数的，不像很多后面的文献中使用上图所表达的带有 $\pi$ 的前缀，当然，是否带有前缀对于模拟方法跟模拟结果并没有什么区别，只是对应的球谐系数会有一些变化。

SH系数的计算就是将输入光源与对应阶数的球谐函数相乘之后进行球面积分，而用离散方式计算则是将环境贴图与球谐函数相乘之后的累加（理论上通过这种方法得到的是env map的球谐系数，也就是光辐射率L的系数，不过由于我们后面使用的是前3阶的SH系数，也就是说高频数据都被滤掉了，按照前面E=AL的说法，可以认为保留下来的数据则只剩下低频的光辐射照度E了，从而可以将 $E_{lm}$ 与 $L_{lm}$ 等同起来）。

这里给出传统积分方法（2小时）与球谐计算方法（1秒钟）的结果对比。

得到光辐射照度E的球谐系数之后，我们就可以通过二次多项式来给出对应位置的光辐射照度E（我们最终渲染需要的是射向视角反向的光照强度，在Lambertian Surface下，这个强度跟E是成正比的）

这个公式是根据前面的球谐函数在笛卡尔坐标系下的格式得到的,具体推导过程是怎么样的呢？从后面给出的数值来看 $c_4=\frac{1}{2\sqrt{\pi}} * \pi = \frac{\sqrt{\pi}}{2},c_2=\frac{\sqrt{3}}{2\sqrt{\pi}} * \frac{2 \pi} {3} /2 = \frac{\sqrt{\pi}}{2\sqrt{3}}$ ，看起来跟后面常用的基函数格式存在一定的关系，相当于 $n_{lm} * c_{lm}/A_l = Coe_{lm}$ ，其中 $Coe_{lm}$ 指的是下图中的各个基函数的前缀系数，比如 $Coe_{00} = \frac{1}{2\sqrt{\pi}}$ ， $A_l$ 则是前面说到过的用于实现从L到E转变的低频滤波函数，其中 $A_{0, 1, 2} = \pi, \frac{2\pi} {3}, \frac{\pi} {4}$ ， $n_{lm}$ 则是上面多项式公式中放在 $c_{lm}$ 之前的系数， $c_{lm}$ 则是上面公式中 $L_{lm}$ 前的c。

公式为：
$E(n) = n^T M n$

$M = \left \{ \begin{matrix} c_1L_{22} & c_1L_{2-2} & c_1L_{21} & 2c_2L_{11} \\ c_1L_{2-2} & -c_1L_{22} & c_1L_{2-1} & 2c_2L_{1-1} \\ c_1L_{21} & c_1L_{2-1} & 3c_5L_{20} & 2c_2L_{10} \\ c_2L_{11} & c_2L_{1-1} & c_2L_{10} & c_4L_{00} - c_5L_{20}\\ \end{matrix} \right \}$

上面公式中 $n$ 是表面法线（列向量） $(x, y, z, 1)^T$ 而 $n^T$ 则是法线的转置（行向量）。

从这里我们可以看到，E(n)的计算可以不再需要Irradiance Map的介入，直接通过代码来实现，计算代码给出如上，只是不知道这里的矩阵M是逐顶点存储，还是整个场景共用，从给出的代码看到矩阵是整个场景共用的，不过RGB三个分量各需要一个矩阵，矩阵跟系数之间的关系给出如下：

    c1 = 0.429043 ; 
    c2 = 0.511664 ; 
    c3 = 0.743125 ; //3*c5
    c4 = 0.886227 ; 
    c5 = 0.247708 ;

    matrix[0][0] = c1*coeffs[8]; // c1 L_{22}  
    matrix[0][1] = c1*coeffs[4]; // c1 L_{2-2} 
    matrix[0][2] = c1*coeffs[7]; // c1 L_{21}  
    matrix[0][3] = c2*coeffs[3]; // c2 L_{11}  

    matrix[1][0] = c1*coeffs[4]; // c1 L_{2-2} 
    matrix[1][1] = -c1*coeffs[8]; //-c1 L_{22}  
    matrix[1][2] = c1*coeffs[5]; // c1 L_{2-1} 
    matrix[1][3] = c2*coeffs[1]; // c2 L_{1-1} 

    matrix[2][0] = c1*coeffs[7]; // c1 L_{21}  
    matrix[2][1] = c1*coeffs[5]; // c1 L_{2-1} 
    matrix[2][2] = c3*coeffs[6]; // c3 L_{20}  
    matrix[2][3] = c2*coeffs[2]; // c2 L_{10}  

    matrix[3][0] = c2*coeffs[3]; // c2 L_{11}  
    matrix[3][1] = c2*coeffs[1]; // c2 L_{1-1} 
    matrix[3][2] = c2*coeffs[2]; // c2 L_{10}  
    matrix[3][3] = c4*coeffs[0]- c5*coeffs[6] ; // c4 L_{00} - c5 L_{20}