首页 > 代码库 > 使用Compute Shader加速Irradiance Environment Map的计算
使用Compute Shader加速Irradiance Environment Map的计算
- Irradiance Environment Map基本原理
Irradiance Environment Map(也叫Irradiance Map或Diffuse Environment Map),属于Image Based Lighting技术中的一种。
Irradiance Map的详细定义可参考GPU Gems 2 Chapter 10.“Real-Time Computation of Dynamic Irradiance Environment Maps”。简单说来就是一种用于近似Environment Diffuse Lighting的方法。想象一个场景中有k个方向光,方向分别为d1…dk,光照强度为i1…ik,对于一个法线和Diffuse Color分别为n和c的Lambert表面,其光照强度为:
对于Environment Lighting,我们可以用一个Cube Map来表示,Cube Map里的每一个texel就是一个方向光,光强度为texel的值,方向为texel的location。这样就能通过一个Cube Map来表示任意的Environment Lighting。一般把这个Cube Map叫做Light Probe。
对于Lambert表面,其光照强度只和法线n和光照方向l相关,所以给定一个Light Probe,可以计算出所有可能的法线方向的光照,然后存储到一个Cube Map里,渲染时,只需要使用法线n去这个Cube Map里索引就能得到Environment Lighting,这个存储着光照的Cube Map就叫Irradiance Map。计算的伪代码如下:
diffuseConvolution(outputEnvironmentMap, inputEnvironmentMap) { for_all {T0: outputEnvironmentMap} sum = 0 N = envMap_direction(T0)
for_all {T1: inputEnvironmentMap} L = envMap_direction(T1) I = inputEnvironmentMap[T1] sum += max(0, dot(L, N)) * I
T0 = sum return outputEnvironmentMap }
对于每一个法线n都需要去遍历所有的光线方向,算法复杂度为O(NM),N为Light Probe的大小,M为Irradiance Map的大小。
- Spherical Harmonics
由于Diffuse光照本身是变化很缓慢的低频数据,所以可以使用SH来加速计算。把算法分为两步:
1. 把Light Probe投影到SH上,求解出SH系数存储下来。
2. 将Light Probe的SH和Diffuse Transfer的SH做卷积即可求出Irradiance Map。
具体做法如下:
//将LightProbe投影到SH上 for each texel of the lightProbe { lightSample = texelRadiance; weight += texelSolidAngle; //计算光照方向 l = texelDirection; //根据光照方向求出SH基函数 SHBasis = calculateSHBasis(l); //累加SH系数 lightSH += lightSample*SHBasis; } lightSH = lightSH*4*PI/weight; for each texel of the irradianceMap { //法线方向 n = texelDirection; // 求出cosine lobe的SH diffuseSH = calculateDiffuseSH(n); // 用cosine lobe的SH和light probe的SH做卷积 irradiance = dotSH(diffuseSH, lightSH); // lambert brdf irradiance *= 1/PI; texelValue = irradiance; }
使用SH来计算的话,Light Probe和Irradiance Map只需要分别遍历一遍,所以算法复杂度为O(KN+KM),N为Light Probe的大小,M为Irradiance Map的大小,其中K为SH系数的个数,对于Diffuse光照,使用3阶的SH函数就能获得不错的近似结果,3阶的SH有9个系数,所以K远小于N和M。
因为Diffuse光照本身是低频的,所以输出的Irradiance Map可以使用较小的分辨率,那么整个算法的开销主要是在第一步——把Light Probe投影到 SH上。
- 使用GPU计算Light Probe SH
GPU Gems 2 的Chapter 10介绍了使用pixel shader来计算Light Probe SH的方法,使用SM5.0的Compute Shader来计算可以获得更大的加速比。
观察求解Light Probe SH的过程——遍历所有的texel,对于每个texel,求解出SH,然后累加,最后累加的结果就是SH系数。如果使用并行的算法,伪代码如下:
g_mutex; for each texel of the lightProbe { lightSample = texelRadiance; //计算光照方向 l = texelDirection; //根据光照方向求出SH基函数 SHBasis = calculateSHBasis(l); //累加SH系数 g_mutex.lock(); lightSH += lightSample*SHBasis; weight += texelSolidAngle; g_mutex.unlock(); } lightSH = lightSH*4*PI/weight;
但是对于GPU的Compute Shader,只能同步一个Group里的thread,对于不同Group的thread无法同步,所以无法使用一个加锁的全局变量不断累加的方法。
既然无法一次求出所有的累加结果,那么就先求出每个Group的累加结果,然后根据GroupID写入到输出的Buffer,然后把这个Buffer作为输入,重复之前的操作,直到输出Buffer的Size为1时就求出了结果。比如对于一个512x512的Cube Map,固定Thread Group大小为8x8,那么我们分配64x64个Group,其输出Buffer大小为64x64(每一个Group输出一个结果),运行Compute Shader计算结果输出到Buffer,这时数据就缩小到了64x64,然后重复之前的操作,下一轮的数据大小就变成了8x8(64/8)。一直重复这个操作,直到输出的Buffer大小为1时就求解出了结果。
这个算法的思路和HDR渲染中求解场景的平均亮度是一样的,在求解平均亮度时,每次把Texture的Size缩小到1/4做Down Sample,直到Texture大小为1x1时就求出了平均亮度。
求解每个Group的SH累加结果的算法参考Nvidia的 “Optimizing Parallel Reduction in CUDA”, Parallel Reduction的思路就是一个递归的tree-based approach,如下图
对于Shader代码,使用循环来模拟这个过程,具体做法是设置一个步长step,把相隔step个步长的数据相加,然后step乘以2,重复这个过程,直到step大于N,N为输入数据的大小,循环累加的代码如下:
for (uint s = 1; s < groupthreads; s *= 2) // stride: 1, 2, 4, 8, 16,32, 64, 128 { int index = 2 * s * GI; if (index < (groupthreads)) sharedMem[index] += sharedMem[index + s]; GroupMemoryBarrierWithGroupSync(); }
算法的流程如下图:
完整的Shader代码如下:
#define THREAD_SIZE_X 8 #define THREAD_SIZE_Y 8 #define GROUP_THREADS THREAD_SIZE_X*THREAD_SIZE_Y Texture2D<float> g_InputBuffer : register(t0); RWTexture2D<float> g_OutputBuffer : register(u0); groupshared float g_ShareMem[GROUP_THREADS]; [numthreads(THREAD_SIZE_X, THREAD_SIZE_Y, 1)] void ReductionCS(uint3 Gid : SV_GroupID, uint3 DTid : SV_DispatchThreadID, uint3 GTid : SV_GroupThreadID, uint GI : SV_GroupIndex) { //加载数据到share memory中 uint Idx = DTid.y*g_InputBuffer.Length.x + DTid.x; g_ShareMem[GI.xy] = g_InputBuffer[Idx]; GroupMemoryBarrierWithGroupSync(); // 循环累加所有数据 [unroll] for (uint s = 1; s < GROUP_THREADS; s *= 2) // stride: 1, 2, 4, 8, 16, 32, 64, 128 { int index = 2 * s * GI; if (index < GROUP_THREADS) g_ShareMem[index] += g_ShareMem[index + s]; GroupMemoryBarrierWithGroupSync(); } if (GI == 0) { //写入结果 g_OutputBuffer[Gid.xy] = g_ShareMem[0]; } }
NV的paper中还提到了可以进一步优化,GPU在运行Thread Group时,会把Thread划分为Warp,Nvidia的GPU中一个Warp包含32个Thread,这些线程是由SIMD32处理器同步运行的,所以如果线程数目小于32时,可以去掉GroupMemoryBarrierWithGroupSync() 的调用来提升性能。(对于AMD的GPU,把线程划分为Wavefronts,和NV的Warp对应,AMD的每个Wavefronts中包含64个同步执行的Thread)
- 运行结果
运行的参考对象是DirectX ToolKit中的SHProjectCubeMap,这个函数使用CPU来计算Light Probe SH。测试使用的Light Probe是一个512x512的R16G16B16A16_FLOAT格式的HDR Cube Map,测试结果发现在Release模式下使用Compute Shader可以比DXTK的CPU版本快10倍以上,加速比很高。
- 其他
1. Light Probe一般都是HDR格式的,所以生成的Irradiance Map也是HDR格式的,那么使用Irradiance Map计算光照时需要使用HDR渲染。
2. Diffuse Lighting一般用3阶的SH足矣,对于HDR的Light Probe,5阶的SH会更准确(4阶的Diffuse Transer SH系数为0,所以无需使用)。当然也意味着更多的计算量。
3. 使用GPU计算的时候,计算的结果参考了DXTK的函数,但是发现和DXTK的结果有些许偏差。仔细调试后发现是浮点累加误差导致的,DXTK的结果并不准确,在其SHProjectCubeMap源码实现中,是使用float类型不断累加,比如对于求解solid angle weight,最后的结果会累加到几百万,但是对于单个的weight数据,其值可能只有1左右,而float有效数字只有6位,所以在累加的过程中会有误差,数据到几百万时,误差可能会达到几十。理论上来讲,solid angle weight是使用texture uv求解的,那么对于Cube Map的每一个face,其累加结果应该都是相同的,DXTK因为使用float累加,浮点误差导致其每个face的solid angle weight累加结果并不相同,所以DXTK的结果并不准确,只能作为参考值。
实际上相比于DXTK,GPU算法的准确度更高,因为调试的时候发现GPU求解出来的结果,每一个face的solid angle累加结果都是相同的,不会像DXTK那样有误差。因为在GPU算法中,累加是分层求解的(参考前面的算法图解),每一层中的节点数据范围都差不多,这样浮点的误差就不会像直接累加那么大,所以GPU求解的速度和精度都好于DXTK的SHProjectCubeMap。
参考资料:
http://http.developer.nvidia.com/GPUGems2/gpugems2_chapter10.html
http://developer.download.nvidia.com/assets/cuda/files/reduction.pdf
http://diaryofagraphicsprogrammer.blogspot.com/2014/03/compute-shader-optimizations-for-amd.html
https://seblagarde.wordpress.com/2012/06/10/amd-cubemapgen-for-physically-based-rendering/
使用Compute Shader加速Irradiance Environment Map的计算