物理光源：Linearly Transformed Cosines

1. 介绍

1️⃣首先，确认LTC技术的目标：解决多边形面光源的光照积分问题，特别是高光部分的求解。多边形面光源的积分和我们之前考虑的问题都不一样：

• 常见光源（点光源、SpotLight等）：简单，积分域就是一个点。
• 环境贴图：复杂一点，积分域是半球域——一般来说，我们都是使用IBL技术类，进行简化。

我们考虑更加物理的光源：多边形光源。这个时候，积分域 $S$ 就很复杂了，求解更加复杂：
$$I={\int }_S(L(w_l)\cdot f(p,w_i,w_o)\cdot \cos ({\theta }_p))d{w}_i$$

注意：我们考虑的是多边形光源，这意味着这个光源的所有顶点都是在一个平面上，不能有凹凸！

最朴素的思想就是：使用蒙特卡罗方法。但这个思路在性能上是无法接受的。

2️⃣作者就发现了一种线性变换球面分布（Linearly Transformed Spherical Distributions ）的思想。这种思想就是对于任意一个球面分布函数，一定可以通过一个线性变换矩阵将其变化到另外一个球面分布函数。这个时候希望来了，由于$co{s}^{+}({\theta }_i)$是一个球面分布函数（余弦分布函数），$f(p,w_i,w_o)$ 也是一个球面分布函数。

那么具体这个算法到底是什么呢？

$co{s}^{+}$是clamp余弦

3️⃣根据论文，我们拿来进行简化计算的球面分布函数，实际上是 $co{s}^{+}({\theta }_i)\cdot \frac{1}{\pi }$。而原来的复杂的球面分布函数，实际上是 $f(p,w_i,w_o)\cdot \cos {\theta }_i$

2. 算法分析

2.1 思路

1️⃣根据第一节，我们有两个认识：

1. 在多边形光源上进行实时光照积分是困难的。($S$ 是多边形积分域)
   $$I={\int }_S(L(w_l)\cdot f(p,w_i,w_o)\cdot \cos ({\theta }_p))d{w}_i$$
2. 根据LTCS思想（数学思想，所以我不懂，只是给出结论）：存在一个变换矩阵 $M$，可以把简单的球面分布函数 $co{s}^{+}({\theta }_i)\cdot \frac{1}{\pi }$，变换成复杂的球面分布函数 $f(p,w_i,w_o)\cdot \cos {\theta }_i$。
   $$f(p,w_i,w_o)\cdot \cos {\theta }_i\approx M\ast \frac{co{s}^{+}({\theta }_i)}{\pi }$$
   注意：这里的*不是普通的矩阵乘法，这里的线性变换是指把我们的入射向量乘以矩阵M。

矩阵$M$ : 三阶矩阵。根据论文，它主要包括如下的性质：

2️⃣但是，还是看不出来有什么用，所以我们需要知道LTCS的另外一个性质：
$${\int }_{S}D(w)dw={\int }_{S_t}{D}_t(w_t)d{w}_t$$
所以，我们会有：
$${\int }_{S}f(p,w_i,w_o)\cdot \cos {\theta }_{i}d{w}_i={\int }_{S_t}\frac{co{s}^{+}({\theta }_i)}{\pi }d{w}_t$$
而积分域：$S=M\ast S_t$。积分域进行变换，实际上就是对这个多边形域的每个顶点，应用这个矩阵$M$。

3️⃣对 $cos$ 进行积分是很简单的，哪怕积分域 $S_t$ 比原积分域 $S$ 复杂很多，也不会有太大影响——实际上，积分域 $S_t$ 也不会比 $S$ 复杂，毕竟，我们并没有增加删除顶点，而只是对顶点的位置进行变换（不考虑裁剪：$S$ 是五边形，$S_t$ 也会是五边形）。而且似乎对 $cos$ 进行积分有近似方法，这个后续进行分析。

4️⃣所以我们的最终的目标就是：放弃在多边形积分域$S$中，使用复杂的 $D(w)=f(p,w_i,w_o)\cdot \cos {\theta }_i$进行光照积分，而是使用简单的 $D_t(w_t)=\frac{co{s}^{+}({\theta }_i)}{\pi }$进行积分计算。
$$I={\int }_{S_t}(L(w_t)\cdot D_t(w_t))d{w}_t$$

但目前还是有两个问题需要解决：

1. 我们需要确定积分域 $S_t$，所以我们需要知道 $M^{-1}$。
2. 光照积分中可不只有球面分布函数 $D_t(w_t)$，还有入射光 $L(w_t)$！

我们依次来解决！

2.2 预计算矩阵M

1️⃣根据原论文，我们可以知道，变换前后的球面分布函数有如下关系：
$$D(w)=D_t(w_t)\frac{\delta w_t}{\delta w}=D_o(\frac{M^{-1}w}{||M^{-1}w||})\frac{|M^{-1}|}{||M^{-1}w|{|}^3}$$
其中，$D(w)=f(p,w,w_o)\cdot \cos {\theta }_i$，$D_o(w_t)=\frac{\cos {\theta }_t}{\pi }$，$w_t=M\ast w$。$M$ 就是我们需要的转换矩阵。根据原论文，这个矩阵 $M$ 实际上有如下固定形式：

2️⃣所以，我们的问题从求得一个矩阵，变成了求得$(a,b,c,d)$的值，那么，怎么求呢？

和IBL方法类似，我们在BRDF中不考虑菲涅尔项，而且假设入射方向等于法线方向，那么 $f(p,w_i,w_o)\cdot \cos {\theta }_i$ 项就只依赖于视线方向和粗糙度。那么我们就可以计算得到一个张类似于BRDF LUT的二维纹理，其横纵坐标是视线法线夹角和粗糙度。

3️⃣咦？我们似乎还是不会求$M$？根据大佬的博客：对于任意一个$f(p,w_i,w_o)\cdot \cos {\theta }_i$ 一定找的到一个M变换矩阵把他变换到余弦分布，因此其实只需要遍历所有的矩阵M总能找到一个误差足够小的矩阵 $M$，不过穷举法不可能，我们采用一种单纯形法用于寻找 $M$。单纯形法的具体过程有点类似梯度下降的思路，我们随机初始化一个矩阵 $M$，然后经过计算，我们会得到一个“梯度”也就是我们该往哪个方向去修正我们的矩阵 $M$，循环往复，直到 $M$ 的误差在接受范围内。

具体过程来说：

1. 对于LUT的每个坐标，也就是一对具体的视线法线夹角和粗糙度，我们可以随意取一个入射方向$w$，求得$N_0=f(p,w,w_o)\cdot \cos {\theta }_i$。
2. 初始化一个矩阵 $M$，代入公式 $D_o(\frac{M^{-1}w}{||M^{-1}w||})\frac{|M^{-1}|}{||M^{-1}w|{|}^3}$求值，得到$N_1$。
3. 求$|N_0-N_1|$的L3范数，作为误差 $e$ 。
4. 不断更新 $M$ ，知道误差 $e$ 足够小（怎么更新，暂未研究）。

4️⃣我们进行预计算，最终得到一张横纵坐标是视线法线夹角和粗糙度，存储了转换矩阵的逆 $M^{-1}$的LUT。

M矩阵预计算code

2.3 恒定多边形光源

1️⃣最简单的情况就是：多边形光源的辐照度是恒定的，那么这个时候，入射光$L$就可以直接提出来。
$$I=L\cdot {\int }_{S_t}\frac{1}{||s-p|{|}^2}\cdot \frac{\cos {\theta }_t}{\pi }d{w}_t$$
其中，$s$ 是渲染点，$p$是多边形光源上的点（变化后的，位于$S_t$上）。而且，根据Daniel的方法，上诉积分具有解析解：
$$\begin{gathered} {\int }_{S_t}\frac{1}{||s-p|{|}^2}\cdot \frac{\cos {\theta }_t}{\pi }d{w}_t=\frac{1}{2\pi }\sum _{j=1}^e{n}_p\cdot \frac{{\Gamma }_j}{||{\Gamma }_j||}{{\rm Y}}_j=I_D \\ -------------------- \\ {\Gamma }_j=Y_j\times Y_{j+1} \end{gathered}$$
其中，$e$ 是多边形光源的边数，$n_p$ 为渲染点 $p$ 法线，$Y_j$ 是着色点到多边形顶点的向量（这些顶点都是经过变换的，位于 $S_t$），${{\rm Y}}_j$为向量之间的夹角（$Y_j$ 和 $Y_{j+1}$之间的夹角 ）。

$\frac{1}{||s-p|{|}^2}$来自于光线衰减。

2️⃣所以，我们最终只需进行如下计算即可：
$$L(p,w_o)=L\cdot I_D$$

2.4 纹理光源

1️⃣纹理多边形光源的辐照度明显是不均匀的，所以我们不能直接将入射光提取出来，但我们可以利用下面这个常用技巧：
$$\begin{gathered} I={\int }_S(L(w_l)D(w_l))d{w}_l\approx I_D\cdot I_L \\ ------- \\ {I}_D={\int }_{S}D(w_l)d{w}_l \\ ------- \\ {I}_L=\frac{{\int }_S(L(w_l)D(w_l))d{w}_l}{{\int }_{S}D(w_l)d{w}_l} \end{gathered}$$
我们把积分拆成两部分，第一部分$I_D$很简单，就是我们所需的完美形式，可以直接使用Daniel的方法求得解析解：
$$\begin{gathered} {\int }_{S_t}\frac{1}{||s-p|{|}^2}\cdot \frac{\cos {\theta }_t}{\pi }d{w}_t=\frac{1}{2\pi }\sum _{j=1}^e{n}_p\cdot \frac{{\Gamma }_j}{||{\Gamma }_j||}{{\rm Y}}_j=I_D \\ -------------------- \\ {\Gamma }_j=Y_j\times Y_{j+1} \end{gathered}$$
2️⃣第二部分 $I_L$ 怎么办呢？作者的方法是，可以看做纹理进行滤波后的值。在预过滤的步骤中，作者使用高斯滤波器，采取不同的过滤核半径进行过滤，对应不同的LOD。

根据论文，这个过滤还有个需要注意的地方：预过滤纹理的值必须定义在纹理空间的每个地方，甚至在纹理之外——作者在纹理周围引入一个边缘，在这个区域，增加了滤镜的半径，使其与纹理相交。

3️⃣对于 $I_L$ ，我们还有最后一个需要解决的问题，那就是用来读取这个纹理的UV和LOD Level。

• LOD Level，论文里面是通过一个一维函数——纹理平面的平方距离$r^2$与多边形的面积$a$之比：$\frac{r^2}{a}$。此外，最终结果需要进行混合，例如，我们得到的LOD Level是3.4，那么我们需要读取Level 3和Level 4的纹理值，然后进行混合（权重各自是0.6和0.4）。
• UV。我们在使用Daniel的方法进行求解的时候，会得到一个指向面光源的向量，这个向量所指面光源上的点是对光照贡献最大的点，记作Q点。我们需要做到就是得到这个Q点，并算出它在这个纹理平面上对应的UV坐标。

4️⃣以下读取纹理贴图的源码来自Monica大佬的Github：

// vPolygonalLightVertexPos: 多边形光源的顶点
// vLooupVector：Daniel方法进行求解得到的向量
vec3 fecthFilteredLightTexture(vec3 vPolygonalLightVertexPos[4], vec3 vLooupVector)
{
	vec3 V1 = vPolygonalLightVertexPos[1] - vPolygonalLightVertexPos[0];
	vec3 V2 = vPolygonalLightVertexPos[3] - vPolygonalLightVertexPos[0];
	vec3 PlaneOrtho = cross(V1, V2);
	float PlaneAreaSquared = dot(PlaneOrtho, PlaneOrtho);
	
	SRay Ray;
	Ray.m_Origin = vec3(0);
	Ray.m_Dir = vLooupVector;
	vec4 Plane = vec4(PlaneOrtho, -dot(PlaneOrtho, vPolygonalLightVertexPos[0]));
	float Distance2Plane;
	rayPlaneIntersect(Ray, Plane, Distance2Plane);
	
	vec3 P = Distance2Plane * Ray.m_Dir - vPolygonalLightVertexPos[0];

	float Dot_V1_V2 = dot(V1, V2);
	float Inv_Dot_V1_V1 = 1.0 / dot(V1, V1);
	vec3  V2_ = V2 - V1 * Dot_V1_V2 * Inv_Dot_V1_V1;
	vec2  UV;
	UV.y = dot(V2_, P) / dot(V2_, V2_);
	UV.x = dot(V1, P) * Inv_Dot_V1_V1 - Dot_V1_V2 * Inv_Dot_V1_V1 * UV.y;
	UV = vec2(1 - UV.y, 1 - UV.x);
	
	float Distance = abs(Distance2Plane) / pow(PlaneAreaSquared, 0.25);
	
	float Lod = log(2048.0 * Distance) / log(3.0);

	float LodA = floor(Lod);
	float LodB = ceil(Lod);
	float t = Lod - LodA;
	vec3 ColorA = texture(u_FilteredLightTexture, vec3(UV, LodA)).rgb;
	vec3 ColorB = texture(u_FilteredLightTexture, vec3(UV, LodB)).rgb;
	return mix(ColorA, ColorB, t);
}

5️⃣最终，我们得到了此时的结果：
$$I=I_L\cdot I_D$$

2.5 考虑菲涅尔项

不知道对不对

1️⃣之前在计算过程中一直忽略了菲涅耳项，因此我们还需要计算菲涅耳项用于对上面计算的光照结果进行矫正。这里的思路也是类似IBL的方法，所以我们还需要使用第二张LUT，也就是大名鼎鼎的BRDF LUT：

然后，实时运行时，直接使用：

// F 就是 F0, 或者说 Specular Color
vec2 envBRDF  = texture(brdfLUT, vec2(max(dot(N, V), 0.0), roughness)).rg;
vec3 specular = prefilteredColor * (F * envBRDF.x + envBRDF.y);

2️⃣所以说，最终结果是：
$$I_{Spec}=(F_0\ast LUT.x+LUT.y)\cdot I_D\cdot I_L$$

2.6 关于漫反射项

1️⃣以上一切都是关于计算高光的，和漫反射无关。

2️⃣而实际上，漫反射的计算要简单的多，因为它没有lob，是均匀的，所以实际上它根本不需要考虑变换矩阵 $M$，在代码上，则是如下：

vec3 LTC_Evaluate(vec3 N, vec3 V, vec3 P, mat3 Minv, vec3 points[4], bool twoSided);
...
vec3 spec = LTC_Evaluate(N, V, pos, Minv, points, twoSided);
spec *= texture2D(ltc_mag, uv).w;
// 矩阵M是单位矩阵，不需要任何转换
vec3 diff = LTC_Evaluate(N, V, pos, mat3(1), points, twoSided); 
col  = lcol*(scol*spec + dcol*diff);
col /= 2.0*pi;

3. 参考

[1]Monica的小甜甜的博客
[2]官方源码
[3][Real-Time Polygonal-Light Shading with Linearly Transformed Cosines]