全局光照技术：Precomputed Radiance Transfer

1. 介绍

物体在不同光照下的表现不同，PRT（Precomputed Radiance Transfer）是一 个计算物体在不同光照下表现的方法。光线在一个环境中，会经历反射，折射， 散射，甚至还会在物体的内部进行散射。为了模拟具有真实感的渲染结果，传统的 Path Tracing 方法需要考虑来自各个方向的光线、所有可能的传播形式，并且收敛速度极慢。PRT使用一种预计算方法，该方法在离线渲染的 Path Tracing 工具链中预计算 lighting以及 light transport， 并将它们用球谐函数拟合后储存，这样就将时间开销转移到了离线中。最后通过使用这些预计算好的数据，我们可以轻松达到实时渲染严苛的时间要求，同时渲染结果可以呈现出全局光照的效果。

预计算部分是PRT 算法的核心，也是其局限性的根源。因为在预计算 light transfer 时包含了 visibility 以及cos项，这代表着实时渲染使用的这些几何信息已经完全固定了下来。所以 PRT方法存在的限制包括：

• 不能计算随机动态场景的全局光照
• 场景中物体不可变动。

以上实际上来自Games202作业2，下面给出我自己的技术理解

个人觉得，Light Map、Irradiance Map、Precomputed Radiance Transfer应该是一脉相承的，都是使用预计算，来解决渲染质量和性能的权衡问题：

1. Light Map。此技术实际上只是把常数项（$f_d$）提取出来，然后没有采用split sum，而是基于实际的对象，没有做任何近似，直接使用路径追踪等复杂离线方法来计算。此技术不仅绑定了光源（静态）、还绑定了对象（静态）。一般引擎中的烘焙过程就是使用了此技术。此技术的限制极大，物体和光源必须都是静态的，一旦变换其一，就需要重新烘焙，来计算Light Map。实际渲染过程中，由于引擎都是以低分辨率计算Light Map，然后将多个物体、多个表面的Light Map存在一张纹理上，所以需要使用特定的LightUV来读取Light Map，然后再乘上$f_d$（或albedo）。
   $$L_{pre}=\int L\cdot V\cdot dot(n\cdot l)dw$$
   

2. Irradiance Map。此技术详情见我的另外一篇博客。这里简单讲下：为了避免Light Map过于严苛的限制，让运动物体可以收益于预计算。Irradiance Map提出了一个重大假设：那就是环境光源无限远，物体上的每一个点，相对于环境光，都是一致的，并且忽略掉自反射和自阴影。这样的话，我们就只绑定了光源，物体在此光源的影响范围内，可以任意的旋转和移动。在预计算公式上，其实也和Light Map没有区别（去掉了V），只是考虑原理不一样了。实际渲染过程中，以法线作为索引，来读取这张环境贴图，然后乘上 $f_d$ 即可。
   $$L_{pre}=\int L\cdot dot(n\cdot l)dw$$

3. Efficient Irradiance Map。此技术详情见我的另外一篇博客。个人目前对于RTR4以及原论文中的关于此技术起源的解释，还是不能理解。所以这里给出我的解释： Irradiance Map低频的特性（相当于低通滤波器），意味着它的入射光照分布可以用光滑曲线去近似，而依然可以取得不错的渲染结果。那么直接给出答案吧（以下都以三阶SH举例数字）：我们可以将光照分布投影到少数几个（27个float）SH系数上，实际使用时重建，这样可以极大压缩内存需求。而这个时候，为了投影到SH上，我们需要使用把原本的积分方程拆成两部分：光照部分L和余弦部分cos，分别对其进行SH投影。这里产生了一个问题：法线n和$L(w)$不同，它独立于w，相当于第二个参数，我们需要对每个可能的方向（总共 $width\ast height$ 个方向）做SH投影，那么就会产生 $27\ast width\ast height$个SH系数，这相当于好几张纹理了，这根本不省内存！幸运的是，受益于SH旋转不变的特性，可以简化为基于法线的缩放系数 $A_l$，来可以克服这个问题，最终的渲染过程也正如我们的期待：计算得到最终的法线$N(\theta ,\varphi )$后，根据法线的球坐标，得到各个球谐基$Y_{l.m}(\theta ,\varphi )$的值，然后代入重建公式：
   $$E(\theta ,\varphi )=\sum _{l,m}{A}_l{L}_{lm}{Y}_{lm}(\theta ,\varphi )$$

4. Precomputed Radiance Transfer。PRT又是另外一个技术路线了，它不关注$\int (L)dw$部分（可以直接按照Efficient Irradiance Map进行SH投影），而是重新捡起了 $V$项，并且还考虑了自反射，也就是说，我们在这关注的是$\int (V\cdot cos)dw$部分，这部分也被称为传输项。这里只是起个头，具体的还是看后文。
   

2. Diffuse Transfer

正如前文所说，PRT技术的核心在于怎么求解转移方程，来解释物体如何对自身投射阴影和散射光线。而对于Diffuse物体来说，应用PRT技术是简单的。

对于光照 $L_p(s)$（整个半球域的光照），我们直接进行SH投影，得到 $n^2$ 个SH参数 $(L_p{)}_i$ 。然后对应的，我们对传输项进行SH投影，也得到 $n^2$ 个SH参数 $(M_p{)}_i$ 。

正如下文所示，随着考虑的要素变多，传输项变得复杂

SH基的正交性提供了一个有用的性质：即给定任意两个球上的函数 $a$ 和 $b$ ，它们的投影满足 $\int a(s)b(s)ds={\sum }_{i=1}^{n^2}{a}_i{b}_i$。换句话说，对带限函数的乘积进行积分，可将其简化为投影系数的点积和。所以：
$$L_o=\frac{\rho }{\pi }\int (L_p(s)\cdot M_p(s))ds=\frac{\rho }{\pi }\sum _{i=1}^{n^2}(M_p{)}_i(L_p{)}_i$$

怎么求SH参数？这里简单提个醒：
$$f_l^m=\int f(s)y_l^m(s)ds=\sum _{i=1}^{i=n}f(s_i)y_l^m(s_i)\Delta s_i$$

知道high level的想法后，现在留下的问题就是：传输项是什么？通过考虑阴影、自反射等因素，我们会有如下三个可选的传输项。

2.1 unshadowed diffuse transfer

这是最简单的传输方程，我们不要考虑阴影和自反射。我们假设 $p$ 处的法线为 $N_p$，那么传输项实际就是钳制的余弦项：
$$M_{N_p}(s)=\max (N_p\cdot s,0)$$
下面是作业中的代码：

for (int i = 0; i < mesh->getVertexCount(); i++)
{
    const Point3f &v = mesh->getVertexPositions().col(i);
    const Normal3f &n = mesh->getVertexNormals().col(i);
    auto shFunc = [&](double phi, double theta) -> double {
        Eigen::Array3d d = sh::ToVector(phi, theta);
        const auto wi = Vector3f(d.x(), d.y(), d.z());
        double cos_ = wi.dot(n);
        //double cos_ = 0.0;
                       
        return max(cos_, 0);
        
    };
    // SHOrder : 使用的SH阶数
    // shFunc : 要投影的函数，这里是转移函数
    // m_SampleCount : 采样点数
    auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
    for (int j = 0; j < shCoeff->size(); j++)
    {
        m_TransportSHCoeffs.col(i).coeffRef(j) = (*shCoeff)[j];
    }
}

2.2 shadowed diffuse transfer

这里我们就要考虑阴影了，其实就是在余弦项的基础上加上可见性 $V$。也就是说，传输项变成了：
$$M_{N_p+Shadow}(s)=\max (N_p\cdot s,0)\ast V_p(s)$$
可见项怎么计算？实际上，一个很简单的思路就是：在当前渲染点的法线半球，进行采样得到方向 $s$ 的同时，沿着这个方向 $s$，进行击中判定。如果击中了，则 $V=0$；否则 $V=1$。

下面是作业中的代码：

for (int i = 0; i < mesh->getVertexCount(); i++)
{
    const Point3f &v = mesh->getVertexPositions().col(i);
    const Normal3f &n = mesh->getVertexNormals().col(i);
    auto shFunc = [&](double phi, double theta) -> double {
        Eigen::Array3d d = sh::ToVector(phi, theta);
        const auto wi = Vector3f(d.x(), d.y(), d.z());
        double cos_ = wi.dot(n);
        //double cos_ = 0.0;
                       
        Ray3f r(v, wi);
        return max(cos_, 0) * (1.0 - scene->rayIntersect(r));
        
    };
    // SHOrder : 使用的SH阶数
    // shFunc : 要投影的函数，这里是转移函数
    // m_SampleCount : 采样点数
    auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
    for (int j = 0; j < shCoeff->size(); j++)
    {
        m_TransportSHCoeffs.col(i).coeffRef(j) = (*shCoeff)[j];
    }
}

2.3 interreflected diffuse transfer

这里进一步考虑了自反射，在考虑传输项之前，我们还需要重写光照方程：
$$\begin{gathered} L_o=\frac{{\rho }_p}{\pi }\int (L_p(s)\cdot M_{N_p+Shadow}(s))ds \\ {L}_o^1=L_o+\frac{{\rho }_p}{\pi }\int \overline{L_p^0}(s)M_{N_p}(1-V_p(s))ds \\ {L}_o^2=L_o^1+\frac{{\rho }_p}{\pi }\int \overline{L_p^1}(s)M_{N_p}(1-V_p(s))ds \\ \cdots \end{gathered}$$
实际上，上诉过程就是物体本身多次互bounce的过程，我们迭代计算多少次，就相当于bounce越多。其中，${\overline{L}}_p(s)$是来自物体 $O$ 自身的、方向 $s$ 的辐射度（radiance）。困难在于，除非入射光来自无限远的光源，否则我们实际上不知道 ${\overline{L}}_p(s)$。

但真的困难吗？$L_o$可以直接计算（虽然这里说是直接，但请注意，我们通篇都是讲的GI，算的是间接光照），而自反射导致的迭代依赖积分似乎很复杂？我们如果把这个过程拆成多个Pass，迭代计算自反射不就成了，反正我们不需要过多考虑性能——这些都是离线计算的。但我们必须注意到一个问题：不管是公式上直接看，还是理论理解，我们在计算传输方程的时候，都需要考虑 BRDF了（对于diffuse，也就是 $\frac{\rho }{\pi }$）！

其实之前的计算也可以考虑 $\frac{\rho }{\pi }$，反正是个常数。但考虑自反射，情况就发生了变化，因为我们的转移函数需要考虑来自物体其他点的辐射度了！此时，考虑 $\frac{\rho }{\pi }$ 是必不可少的！

$$\begin{gathered} M_{N_p}(s)=\frac{{\rho }_p}{\pi }\max (N_p\cdot s,0) \\ ------------------- \\ {M}_{N_p+Shadow}(s)=\frac{{\rho }_p}{\pi }\max (N_p\cdot s,0)\ast V_p(s) \end{gathered}$$
以上就是自反射情况下，对之前两个传输方程的修正。那么，让我们正式开始描述把！

• Pass 1。按照shadowed diffuse transfer描述的那样，计算得到模型上每一个点的传输参数——$(M_p{)}_i^0$ 。
• Pass 2。也是和Pass 1一样，遍历模型的每一个顶点。以 $p$ 点为例，在其法线半球域，进行采样。对于采样得到的方向 $s_d$，我们只考虑和物体本身相交的，也就是说，传输方程参数的计算如下：
  $$(M_p{)}_i^1=\sum _{s_d}(1-V_p(s_d))\cdot (M_q{)}_i^0\cdot M_{N_p}(s_d)$$
  其中，$(M_p{)}_i^1$ 的下标代表球谐参数的索引（第几个SH参数），其上标代表迭代次数。$q$ 代表物体上和 $s_d$相交的点。上述方程比较好理解：$1-V_p(s_d)$ 是由于我们只考虑和物体相交的方向；$(M_q{)}_i^0$则是需要考虑 $q$ 点本身的光线传输情况。

$(M_q{)}_i^0$ 实际上并不是 $q$ 的传输参数，或者说，$q$ 点身可能并不是模型的一个顶点，而是相交的三角形图元上的一个点。所以$(M_q{)}_i^0$ 实际上是这个相交图元的三个顶点的传输参数的插值。Games202里面是说的是重心坐标插值。

• Pass 3，就是迭代，直接贴公式了：
  $$(M_p{)}_i^2=\sum _{s_d}\frac{{\rho }_p}{\pi }\cdot (1-V_p(s_d))\cdot (M_q{)}_i^1\cdot M_{N_p}(s_d)$$
• 可能继续的Pass4、Pass 5。只要你不嫌麻烦。
• Pass final，组合直接传输和多次传输（再说一遍：我们通篇都是讲的GI，算的是间接光照）：
  $$(M_p{)}_i=(M_p{)}_i^0+(M_p{)}_i^1+(M_p{)}_i^2+\cdots +(M_p{)}_i^n$$

这里就暂时不贴代码了，嘛哈哈哈，当时太菜了，写错了。

2.4 实时计算

这个就很简单了，以三阶SH投影为例：

• light。对于光照的投影（例如：来自环境贴图），我们得到 $9_R+9_G+9_B=27$ 个SH参数，然后实际上可以以uniform mat3传入顶点着色器。由于是27个参数，我们应该传入 $Ma{t}_R,Ma{t}_G,Ma{t}_B$ 三个mat3矩阵。
• Transform。对于转移函数的投影，则是有 $Numbe{r}_{顶点数}\ast 9$ 个SH参数。所以我们应该把这个转移SH参数，设置成顶点数据结构体的成员，类似于顶点法线、顶点位置等。只不过，这个成员是mat3类型。

2.5 旋转光源

如果我们想要旋转光源，那就破坏了预计算系方法的静态光源要求，除非我们重新预计算旋转后光源的分布情况。但更好的方法还是在SH系数上花功夫。本节介绍这个方法，它主要利用了球谐函数的两个性质：

• 球谐函数具有旋转不变性，通俗的讲就是：假设有一个旋转 $R$，对原函数 $f(x)$ 的旋转 $R(f(x))$ 与直接旋转$f(x)$ 的自变量是一样的，即 $R_1(f(x))=f(R_2(x))$。

• 对每层band上的SH coefficient，可以分别在上面进行旋转，并且这个旋转是线性变化的。 也就是说，SH系数可以看成是向量，并且可以拆分。

⭐️ 我们对环境贴图（举个例）进行SH投影，得到9个SH系数，可以把这个看出一个$vec9$向量，旋转环境贴图实际上就是旋转这个vec9向量：
$$\begin{gathered} Light\_SH\_vec9=[L_{0,0},\cdots ,L_{2,2}]=[\sum L_{n_p}\cdot Y_{0,0}(n_p)\cdot \Delta Area,\cdots ,\sum L_{n_p}\cdot Y_{2,2}(n_p)\cdot \Delta Area] \\ ----------------------- \\ R(Light\_SH\_vec9)=R([L_{0,0},\cdots ,L_{2,2}])=R(\sum L_{n_p}\cdot \Delta Area\cdot [Y_{0,0}(n_p),\cdots ,Y_{2,2}(n_p)]) \\ =\sum L_{R^{-}(n_p)}\cdot \Delta Area\cdot R([Y_{0,0}(n_p),\cdots ,Y_{2,2}(n_p)]) \end{gathered}$$
通过如上公式推导，我们可以知道：$R(Light\_SH\_vec9)$的求解可以看作 $R([Y_{0,0}(n_p),\cdots ,Y_{2,2}(n_p)])$的求解（至于$R^{+}(L_{n_p})=L_{R^{-}(n_p)}$，就是换个采样的反射度值，是一个未知且无关的变换）。而$R(\cdot )$实际上就是一个$mat9$矩阵（拿来旋转 $vec9$ 啊！）。此外，$[Y_{0,0}(n_p),\cdots ,Y_{2,2}(n_p)]$ 实际上就是法线$n_p$的SH投影，或者说 $P(n_p)=[Y_{0,0}(n_p),\cdots ,Y_{2,2}(n_p)]$。最终，我们需要考虑的就是$R(P(n_p))$的求解!
$$R_1(P(n_p))=P(R_2(n_p))$$

1️⃣首先，以三阶SH为例，来确定符号。根据上面的分析，投影函数$P(n_i)$ 是一个类球谐函数（多个球谐函数组成的向量，从下面的代码也可以看出来），返回SH系数向量，返回类型暂定为 $vec9$。它的伪代码大致如下：

/*float getSH(int l, int m)
{
	float sh;
	for(int i = 0; i < width; ++i)
	{
		for(int j = 0; j < height; ++j)
		{
			vec3 n =GetNormal(i, j)
			sh += envMap(i, j) * Y(l, m, n) * getArea(i, j);
		}
	}
	return sh;
}*/

vec9 P(vec3 n)
{	
	vec9 sh;
	for(int l = 0; l < 3; ++l)
	{
		for(int m = -l; m <= l; ++m )
		{
			sh[getID(l, m)] = Y(l, m, n);
		}
	} 
}

所以 $P(n_i)$具备旋转不变性，所以有如下关系：
$$M_{mat9}\cdot P(n_i)=P(R(n_i))$$
我们知道$P(n_i)$ 以及 $P(R(n_i))$的值（后者就是对法线进行旋转），似乎不好求解得到$M_{mat9}$。但如果我们采样9个法线，并把他们组合成矩阵 $A=[P(n_0),P(n_1),\cdots ,P(n_9)]$，根据数字知识，此时满足：
$$M_{mat9}\cdot A=[P(R(n_0)),P(R(n_1)),\cdots ,P(R(n_9))]$$
这个时候就好办了，A是矩阵，我们就可以利用矩阵的逆，来进行操作了：
$$M_{mat9}=[P(R(n_0)),P(R(n_1)),\cdots ,P(R(n_9))]\cdot A^{-1}$$

2️⃣接下来可以讨论具体算法的实现了：

1. 首先，选取9个法线 $[n_0,n_1,\cdots ,n_9]$，对其中的法线$n_i$，在球谐上投影，得到$P(n_i)$。9`个 $vec9$ 向量$P(n_i)$构成了我们需要的矩阵 $A$，并求出 $A^{-1}$。
2. 其次，给定旋转$R$，对所有 $n_i$ 依次做旋转 $R(n_i)$、SH投影取系数$P(R(n_i))$。 此时我们应该得到9个 $vec9$ 向量 $P(R(n_i))$，这些向量构成了我们需要的矩阵$S$。
3. 接着，求出所需的矩阵 $M=S{A}^{-1}$，用 $M$ 乘以Light SH coefficient vector就可以得到旋转后的SH系数了。

3️⃣Games202里面给出了优化思路：mat9太大了，且计算耗时，而且0阶那个系数根本不需要考虑旋转。在此基础上，我们可以把剩余的8个系数，按照所属阶拆成两部分：1阶的vec3，2阶的vec5。具体流程和上面的类似，可以看看如下截图：

4️⃣不管实现方法如何，我们最终可以发现，这个用来旋转光源SH Vector的函数（一个，或多个矩阵），实际上是独立于环境光的。也就是说，我们完全可以预计算几个常用的旋转角度，通过上诉方法求出旋转矩阵，然后实际运行的时候就可以直接用了。其次，这个旋转算法不是是适用于这个PRT，它依赖的仅仅是球谐函数的特性，所以其他使用SH基进行压缩的GI算法，也可以使用这个方法来旋转光源。

3. Specular Transfer

对于Specular物体来说，应用PRT技术则有点复杂。我目前也没有完全弄明白，这里先给出我们需要提前知道的设定：

• BRDF方程不依赖法线，而依赖于反射向量$R$：$f_r=G(s,R,r)=BRDF(s,R,roughness)$。
• BRDF不是各向异性的，所以它的lobe是以$R$为轴，圆对称的。

3.1 shadowed specular transfer

我们首先给出渲染方程：
$$L_{GS}(L_p,R,r)=\int L_p(s)G(s,R,r)V_p(s)ds$$
如果依然按照之前的思路，会遇到困难：$G$不仅是入射光方向$s$的函数，不能被简化成SH系数向量。

作者的思路是：

1. 首先，去除$R$的影响，让：$G_r^{\ast }(z)=G(s,(0,0,1),r)$
2. 然后，将积分看作是：$G_r^{\ast }$和光照进行卷积，然后投影到SH基上，最后代入实际的$R$，来获得$L_{GS}(L_p,R,r)$。
3. 重建公式：$L_{GS}=\sum {\alpha }_l^0\cdot G_l^0\cdot L_l^m\cdot Y_l^m(R)$

SH Convolution一个圆对称核函数 $h(z)$，与函数 $f$ 的卷积表示为：$h\ast f$。卷积的投影满足：
$$(h\ast f{)}_l^m=\sqrt{\frac{4\pi }{2l+1}}{h}_l^0{f}_l^m=a_l^0\cdot h_l^0\cdot f_l^m$$

1️⃣综上所诉，我们首先可以把 $G(s,(0,0,1),r)$ 投影到SH基上。方法如下：
$$G_l^m=\int G(s,(0,0,1),r)y_l^m(s)\mathrm{d}s$$

2️⃣

ToDo