如何利用LTC实现实时多边形面积光？

2021-03-18

2814

要点概述及概念理解

论文的目的是要实现实时的多边形面光源的反射，并且能有基于物理的brdf效果。而要实现基于物理的brdf反射，我们要计算每个着色点来自面光源的irradiance。这里有两个难点：

1.物理brdf球面的多边形积分没有解析解，而且非常消耗。

2.面光反射需要带有真实感物理的效果。

为了解决这两个问题，我们引入一种新的球面分布叫linearly transformed spherical distribution(LTSDs)。简单来说就是一个球面分布用一个3x3的矩阵做了线性变换后得到一个新的分布。经过这个线性变换后，改变了原分布的“形状”，例如roughness，elliptic anisotropy，skewness。我把论文的图片贴出来，但必须让大家明白，该图片描述的是LTSDs，并不是后面提到的LTC(Linearly Transformed Cosine)。

下面引入面光反射用到的分布，叫做Linearly Transformed Cosine(LTC)，注意，这个分布是LTSDs这类分布中的其中一员而已。下面的图是4种LTSDs，其中包括第三列的LTC。

经过线性变换后的分布，继承了原分布的一些特点，例如归一化，球面多边形积分，重要性采样等。

而这些分布，都可以改变形状近似出基于物理的brdf，为何我们要选择原分布是clamped cosine呢？

核心思想

介绍完基本的概念后，我整理了整个实现的核心思想。

1.在实时渲染的每个pixel shading中，我们利用LTC这个分布，去做面光的反射，也就是说我们要在shading中计算出着色点的irradiance。

2.由于LTC的分布可以转换回原来的cosine分布，在原分布中可以求得多边形的irradiance的积分解析解，这个irradiance等同于在LTC分布中的面光积分。

3.这个多边形在半球上的irradiance的积分的算法，详细参考Geometric Derivation of the Irradiance of Polygonal Lights。

4.上面提到了从原cosine分布转换到LTC有一个3x3矩阵M，这个矩阵M需要通过一个离线算法去生成，可以达到LTC拟合出GGX分布等效果，后面详细介绍。

5.各个分布之间的关系：

这里可能会产生一个疑问，θv应该是出射光才对，为何论文说是入射光？由于BRDF是双向的，也就是说，出射光和入射光交换可以得到同一个效果，我在实时shading的时候，入射光作为出射光的方向来查询能得到相同的效果。

到现在我们可以知道，我们需要的是两个参数α和θv，查询abcd四个值，所以查询表可以用一张2D纹理来存储数据。

然而，如何计算出矩阵M的abcd四个参数，论文并没提及，幸好作者开源，详细的预计算过程只能自己阅读源码。

下面是我根据源码，总结出的计算思路。

1.估算brdf在出射光方向为θ，粗糙度为α时的reflectance和对应的入射光方向，reflectance这个概念参考pbr-book 8.1节，采样的入射光的反射率除以reflectance做一个概率密度用。

2.计算LTC拟合brdf的矩阵M（，由于LTC是近似在θ和α下的brdf的reflectance，那么这个近似必然存在误差值，为了使误差值最小，源码作者使用了一种叫downhill simplex method的方法去计算矩阵M的四个参数。

3.downhill simplex method算法需要做第一次猜测，然后迭代求使得误差值最小。第一次猜测用的是brdf估计出来的入射光方向做矩阵旋转，每一次迭代用当前的矩阵M，通过重要性采样brdf和LTC做误差，利用蒙特卡洛积分求出总的误差值。

4.上面提到了重要性采样求误差值，在源码中实际就是要采样某个入射光的贡献值，而LTC的采样就利用到公式(1)去求得重要性采样出来的入射方向L的反射率，再和近似的BRDF的L方向的反射率做一个差值进行误差比较。

所以论文说了很多但并没有把公式(1)的应用写出来，而是在作者的源码中得到了体现。

实时着色

理解清楚了LTC的特点后，实时着色就相对比较简单了。我们已经离线生成了构建矩阵M的纹理，那么根据View向量的Zenith Angle和着色表面的粗糙度，查询出对应的M矩阵的参数，并构建出M的逆矩阵。

M的逆矩阵是为了让LTC变换回Do的cosine分布做多边形的积分求出irradiance。

下面是矩形面积光的irradiance计算，请留意哪个orthonormal的坐标系，这里我有个问题想不明白，请看代码注释：

// construct orthonormal basis around N
half3 T1,
T2;
T1 = normalize(V - N * dot(V, N));
T2 = cross(N, T1);

// rotate area light in (T1, T2, N) basis
// 矩阵要变换到T1 T2 N的坐标系，应该是以它们作为基的逆矩阵，这里为何是要乘以原矩阵呢？
// 希望有高人能解答一下，源码是transpose的也就是逆矩阵。
Minv = mul(Minv, float3x3(T1, T2, N));

half3 L[5];
L[0] = mul(Minv, points[0].xyz - P);
L[1] = mul(Minv, points[1].xyz - P);
L[2] = mul(Minv, points[2].xyz - P);
L[3] = mul(Minv, points[3].xyz - P);
L[4] = L[0];

int n = 4;
//下面这个是计算多少个顶点在平面下面，一般不需要做
//ClipQuadToHorizon(L, n);
if (n == 0) return half3(0, 0, 0);

// project onto sphere
L[0] = normalize(L[0]);
L[1] = normalize(L[1]);
L[2] = normalize(L[2]);
L[3] = normalize(L[3]);
L[4] = normalize(L[4]);

// integrate
float sum = 0.0;

sum = IntegrateEdge(L[0], L[1]);
sum = IntegrateEdge(L[1], L[2]);
sum = IntegrateEdge(L[2], L[3]);
if (n >= 4) sum = IntegrateEdge(L[3], L[4]);
if (n == 5) sum = IntegrateEdge(L[4], L[0]);

sum = twoSided ? abs(sum) : max(0.0, sum);

half3 Lo_i = half3(sum, sum, sum);

至于多边形在半球上的irradiance积分，可以参考文章Geometric Derivation of the Irradiance of Polygonal Lights

(https://hal.archives-ouvertes.fr/hal-01458129)。

根据这个算法，矩形的顶点顺序需要逆时针或顺时针传递到shader里。

oid UpdateRectPoints() {
float halfX = transform.localScale.x * 0.5f;
float halfY = transform.localScale.y * 0.5f;
rectPointsInWorld[0] = transform.position - transform.right * halfX - transform.up * halfY;

rectPointsInWorld[1] = transform.position transform.right * halfX - transform.up * halfY;

rectPointsInWorld[2] = transform.position transform.right * halfX transform.up * halfY;

rectPointsInWorld[3] = transform.position - transform.right * halfX transform.up * halfY;

Shader.SetGlobalVectorArray("_RectPoints", rectPointsInWorld);

最后的效果图：