计算机图形学：使用体素锥体跟踪实现全局渲染
似乎都需要的头图
细节没怎么扣，有些摩尔纹
前言
每次面试的时候都需要介绍一下自己的工作，有可能是自己的表达能力的关系吧，面试官似乎并不太在意我做的东西，有时候听得云里雾里的。所以自己写一下整个工作的具体实现，既是帮助自己锻炼语言能力，也是求关注，希望得到潜在的hr认可，所以写下global illumination using voxel cone tracing 的全部实现细节。所有的实现代码都在gitee上公开，写得很烂。。比较关注实现，没怎么进行过优化
gitee:
https://gitee.com/zayrq/svo_renderer/tree/master
视频：
https://www.bilibili.com/video/BV1pa411E7nv?spm_id_from=333.999.0.0&vd_source=9ee6e6b6a85ae4996c878c6bec6ec0cb
技术总览
这是对论文Interactive Indirect Illumination Using Voxel Cone Tracing的一个实现。使用稀疏八叉树对场景（sponza）进行建模，把场景的信息进行体素化（voxelization），并把每个节点的体素值（albedo），保存在一个巨大的3D纹理中。并对纹理的节点进行自底向上的mipmapping。利用四线性插值和mipmapping，实现对cone tracing的近似。使用稀疏八叉树（sparse octree）来加速gpu的计算。
前面的介绍，看起来似乎是一个挺复杂的工作。原理其实非常简单，其实就是把渲染场景的信息保存在一个很大的3D纹理中，然后对3D纹理进行mipmapping(glGenerateMipmap())，实际渲染的时候，核心算法就是一个ray marching的过程。对图形学比较熟悉的朋友应该很容易理解。
// 核心的cone tracing过程
// ori 坐标
// dir tracing 的方向
// tmin cone tracing的最小距离
// tmax cone tracing的最远距离
// near 锥体的高
// aperture 锥体的半径
// VOXEL_SIZE = 1.0 / 512.0
// INV_VOXEL_SIZE = 512.0
// EPS = 1.0 / 128.0
vec4 coneTracing(float tmin, float tmax, vec3 ori, vec3 dir, float near, float aperture)
{
float aperNear = (2.0 * aperture) / near;
vec4 alpha = vec4(0.0);
float width = aperNear * tmin;
for (; tmin < tmax && abs(alpha.a - 1.0) > EPS; tmin += width * 0.5)
{
width = aperNear * tmin;
vec3 pos = ori + tmin * dir;
if (IsIllegalPos(pos) == 1)
break;
float lod = clamp(log2(width * INV_VOXEL_SIZE), 0.0, 7.0);
vec4 a = textureLod(voxelTexture, pos, lod);
alpha += (1.0 - alpha.a) * a
}
return alpha;
} 这里的关键是。可能绝大部分的像素都是无效像素(0.0)，只有少部分存在建模信息的地方是有效的。如果把所有的信息都保存在一个3D纹理中的话，会导致纹理变得无比巨大，纹理还需要生成mipmapping。一个512x512x512的GL_RGBA8纹理，需要512MB的空间保存，这还不包括mipmapping的空间，如果采用更细粒度来保存场景信息，需要的空间更大，这里边绝大部分都是无效的像素。我观察过git上，有不少的实现都是采用这个方法的(git上还有个是用clipmap实现的改良版，不过就不介绍了)。所以，论文提出使用稀疏八叉树来对场景进行体素化，把所有没有信息的节点都直接去掉，只保存有建模信息的节点。对于sponza场景来说，brickSize=3，resolution=512。最终生成的3D纹理，大小不过70MB，并且包括所有的mipmapping结果。远远比之前提到的做法要优秀。但是，使用系数八叉树并不是没有问题。使用稀疏八叉树生成的纹理，像素的坐标无法反映场景的实际位置，必须通过遍历八叉树来查找到对应的像素。而且，必须手动对3D纹理进行mipmapping，不能依赖opengl提供的box filter来一键实现。在tracing的过程中，也需要自己实现textureLod函数，这里面存在大量坑，在网上很难搜索到对应的解决方法（可能是我太菜了吧）。后续的内容将会围绕场景体素化，八叉树构建，isotropic mipmapping, anisotropic mipmapping，cone tracing的具体实现展开。
技术细节
场景体素化(voxelization)，保守化渲染
要把场景的信息保存到3D纹理中，首先的首先，是把整个场景的三角片面，都光栅化成片元，利用片元中的值（位置，颜色），来生成稀疏八叉树和3D纹理。所谓的体素化，实际就是利用显卡的渲染管线，对场景进行光栅化操作，并把光栅化生成的片元，保存在一个SSBO中，而不是输出到纹理。光栅化需要尽量保证生成的片元足够多，可以覆盖整个三角形。这里需要在几何着色器做两件事情：主轴投影，以及保守化渲染(conservative rasterization)。主轴投影需要在几何着色器阶段，把三角形投影到主轴上，所谓的主轴，指的是三角形法线的三个分量中，绝对值最大的那个轴：
uint getMainAxis(vec3 norm)
{
vec3 absNorm = vec3(abs(norm.x), abs(norm.y), abs(norm.z));
if (absNorm.x > absNorm.y && absNorm.x > absNorm.z)
return 0;
else if (absNorm.y > absNorm.z)
return 1;
return 2;
}
vec2 projectToAxis(vec3 pos, uint axis)
{
if (axis == 0)
return pos.yz;
if (axis == 1)
return pos.xz;
return pos.xy;
} 主轴投影可以保证投影的三角形，面积最大。而保守化渲染策略则是基于渲染管线一个渲染逻辑：在光栅化阶段时，只有像素的中心被三角形覆盖，才会生成对应的片元。如果三角形过小，无法覆盖到像素的中心时，此时不会生成片元，会导致最终生成的体素不连续，所以，为了保证三角形覆盖到的像素都能生成片元，需要把三角形的三条边向法线方向扩张：
// 计算边的法线
vec3 GetLine(vec2 a, vec2 b, vec2 center)
{
vec2 diff = normalize(a - b);
vec2 cb = center - b;
vec3 ret;
ret.xy = normalize(dot(cb, diff) * diff + b - center);
ret.z = -(ret.x * b.x + ret.y * b.y);
return ret;
}
// 重新计算三角形的顶点
vec2 GetIntersection(vec3 lineA, vec3 lineB)
{
vec2 ret;
ret.y = -(lineB.x * lineA.z - lineA.x * lineB.z) / (lineB.x * lineA.y - lineA.x * lineB.y);
ret.x = -(lineA.y * ret.y + lineA.z) / (sign(lineA.x) * (abs(lineA.x) + 0.0001));
return ret;
}
// 主轴投影
vec3 pos[3];
pos[0] = vec3(projectToAxis(gl_in[0].gl_Position.xyz, mainAxis), gl_in[0].gl_Position.w);
pos[1] = vec3(projectToAxis(gl_in[1].gl_Position.xyz, mainAxis), gl_in[1].gl_Position.w);
pos[2] = vec3(projectToAxis(gl_in[2].gl_Position.xyz, mainAxis), gl_in[2].gl_Position.w);
// 找到三角形的中间点
vec3 center = (pos[0] + pos[1] + pos[2]) / 3.0;
vec3 lines[3];
// 计算三角形的三条边的法线
lines[0] = GetLine(pos[1].xy, pos[0].xy, center.xy);
lines[1] = GetLine(pos[2].xy, pos[1].xy, center.xy);
lines[2] = GetLine(pos[0].xy, pos[2].xy, center.xy);
// 法线向外扩张
lines[0].z += VOXEL_SIZE * 2.0;
lines[1].z += VOXEL_SIZE * 2.0;
lines[2].z += VOXEL_SIZE * 2.0;
// 重新计算三角形的顶点
pos[0].xy = GetIntersection(lines[2], lines[0]);
pos[1].xy = GetIntersection(lines[0], lines[1]);
pos[2].xy = GetIntersection(lines[1], lines[2]); 
// |-----------|-----------|-------| // 32----z-----20----y-----8-x-0-8-0 // |---|-------|-------|-------|---| // 32--|---b---|---g---|---r---|bz-| // |---|---|-------|-------|-------| // |by-|bx-|--nz---|--ny---|--nx---| // 代码节选 // 获取纹理的颜色 vec4 texColor = texHandles[fs_in.texInd] != invalidHandle ? texture(sampler2D(texHandles[fs_in.texInd]), fs_in.texcoord).rgba : vec4(fs_in.color, 1.0); uvec3 color = uvec3(texColor.rgb * 256.); color = clamp(color, uvec3(0u), uvec3(255u)); // 对片元的坐标进行编码 uvec3 voxelPos = clamp(uvec3(fragPos), uvec3(0u), uvec3(voxelResolution - 1.0)); fragmentlist[cur].x = ((voxelPos.x >> 4u) & 0x000000ffu) | ((voxelPos.y & 0x00000fffu) << 8u) | ((voxelPos.z & 0x00000fffu) << 20u); fragmentlist[cur].y = ((voxelPos.x & 0x0000000fu) << 28) | (color. x << 4u) | (color.y << 12u) | (color.z << 20u); // 对法线编码 vec3 normal = normalize(fs_in.normal); normal = normalize(fs_in.normal); fragmentlist[cur].z |= encodeVec3ToUint(normal);
对片元的处理也非常简单，只是把片元的坐标，法线，片元的颜色值保存在一个SSBO中而已。
八叉树构建
前面我们已经获得整个渲染场景生成的片元了，那么在拥有这些片元之后，利用片元中编码的坐标信息，生成稀疏八叉树就变得非常简单了。只需要遍历所有的片元，利用片元中的坐标信息，自顶向下地生成整个八叉树。我们先把八叉树节点的数据结构讨论一下。
uvec2 octrees[] // |--|------------------------------| -- x // |ab|-----------childInd-----------| // |--|------------------------------| -- y // |cd|-----------brickInd-----------|
八叉树的节点非常简单，每个节点是一个uvec2的向量(简称为node)。其中a、b、c、d都是标记位，childInd是子节点块的地址，brickInd是节点对应的brick在3D纹理中的地址。注意，childInd是子节点块的地址，一个块里面有8个子节点，也就是八叉树的所有八个子节点是连续保存的。brick的设定在下一章会进行解释，当前生成八叉树的阶段，还不需要用到这个数据。这里再解释标记位，本次生成八叉树的过程会使用到的标记位是a、b。a表示节点是否为中间节点，b表示节点是否为叶节点（叶节点的childInd为0，没有子节点）。在稀疏八叉树生成的阶段，我们会通过标记节点的a、b两个标记位，判断是否需要生成当前节点，以及它的子节点。讲到这里，应该就很好理解该如何生成八叉树了吧。按照自己需要的八叉树层级（我使用了8层），对每个片元，都进行一次八叉树遍历，如果节点没有子节点（当前还没生成子节点地址），则对a/b进行标记，并终止遍历。然后对所有被标记的节点，生成子节点块地址（childInd），以及一个纹理块地址（brickInd），这通过使用一个原子计数器来完成。这样子循环8次，直到整个八叉树构建完成。
// level = 8
// fragPos为片元中的坐标，因为我构建的八叉树有8层，所以fragPos的值域就是[0, 255]
// curLev = 0
for (int scale = (1 << (level - 1)); curLev < level; ++curLev, scale >>= 1)
{
// 获取片元位于的子节点偏移，等价于
// idx.x = fragPos.x < scale ? 0 : 1;
// idx.y = fragPos.y < scale ? 0 : 1;
// idx.z = fragPos.z < scale ? 0 : 1;
ivec3 idx = ivec3(step(ivec3(scale), fragPos));
// 获得子节点的地址
ind = (octrees[ind].x & 0x3fffffffu) + (idx.x | (idx.y << 1) | (idx.z << 2));
// 减去scale，获得在新节点下的位置，等价于
// fragPos -= scale * idx;
fragPos &= ~scale;
// 当前节点没有子节点，离开循环
if ((octrees[ind].x & 0x3fffffffu) == 0x00000000u)
break;
}
if (curLev == level)
octrees[ind].x |= 0x40000000u;
else
octrees[ind].x |= 0x80000000u; 核心的遍历代码也是非常简单的，用了位运算代替减法运算，用了opengl函数来代替判断，主要是为了提高计算速度。
纹理的isotropic/anisotropic mipmapping
上一章介绍了整个稀疏八叉树的构建过程，这一章，将要介绍的就是SVOGI中最核心的内容，light injecting（光线插入），以及mipmapping的过程。首先的首先，我们需要对场景中的遮挡信息进行建模（alpha），制作一个只有alpha通道的3D纹理。在完成光线插入的工作之后，我们需要把这个alpha Texture的数据复制到lighting Texture的alpha通道中。完成alpha texture的构建之后，再进行光线的插入（light injecting）。实际上，二者的操作是完全一样的；只是alpha Texture在完成之后基本不会改变，而lighting Texture会频繁改变，需要实时完成。
Brick
我们在上一章中讲过的稀疏八叉树的数据结构，其中，每个节点除了保存一个指向子节点块的地址之外，还会保存一个指向3D纹理中的地址索引，并且在上一章，我们也完成了对brickInd的赋值。SVO中的每个节点，都指向一个3x3x3的纹理块，这个被称为Brick。这些Brick，保存着整个场景的信息。所谓的对场景信息的建模，就是对这些Brick中的值进行更新。每个中间节点的Brick，对应着它的子节点的Mipmapping结果。一般的，坐标的值会保存该节点自身的信息。而剩下的值，则保存纹理的边界信息，也就是邻近节点的纹理信息，这个邻近节点的信息非常重要，由于3D纹理中的邻近块，在几何上其实不是连续的，为了加速计算，也是为了利用gpu的线性插值功能，我们需要把节点的边界信息，保存在节点的块上，这样在tracing的时候，就不需要访问邻近节点，也不需要手动编写三线性插值的代码，会极大地加速计算。这个离散保存的3D纹理，以及线性插值的问题，就是SVOGI在实现中最容易忽略的问题，如果不注意，即使实现了Mipmapping，在最后也得不到正确的结果。
light injecting
光线/场景信息插入的操作过程非常简单，核心的代码跟上一章中对八叉树进行标记是一样的，只是这次不再需要对标记位进行标记而已：
struct BrickPool {
layout(bindless_image, r32ui) volatile uimage3D albedos;
layout(bindless_image, rgba8) image3D lights;
layout(bindless_sampler) usampler3D lightMap;
uint begin;
uint num;
uint brickDem;
int depth;
};
// 当初的设计冗余，实际上不需要数组
layout(binding = 5, std430) buffer Pool
{
BrickPool brickPool[];
};
uint GetBrickBlock(uint brickInd)
{
if (brickInd < brickPool[0].num)
return 0u;
return 1u;
}
// 设计冗余，对纹理的大小设计的比较好的话，比如384*384*depth的话，直接对brickInd做位运算就可以得到纹理的坐标
ivec3 GetBrickCoord(uint brickInd, ivec3 texcoord)
{
uint blockInd = GetBrickBlock(brickInd);
uint brickDem = brickPool[blockInd].brickDem;
uint ind = brickInd - brickPool[blockInd].begin;
uint z = ind / (brickDem * brickDem);
uint y = ind % (brickDem * brickDem);
uint x = y % brickDem;
y /= brickDem;
ivec3 ret = ivec3(x, y, z);
ret *= 3;
return ret + texcoord;
}
uvec4 AccumulateColor(uvec4 oldColor, uvec3 color)
{
float invNum = 1.0 / float(oldColor.a + 1);
uvec4 res = uvec4(0u);
oldColor.rgb = oldColor.rgb * oldColor.a;
res = uvec4(oldColor.rgb + color, oldColor.a + 1);
res.rgb = uvec3(vec3(res.rgb) * invNum);
return res;
}
uvec3 ExtractBrickLocFromFragment(uint fragInd)
{
uvec3 res = uvec3(0u);
res.z = fragmentList[fragInd].y & 0x0000000fu;
res.y = fragmentList[fragInd].z >> 28u;
res.x = (fragmentList[fragInd].z & 0x0f000000u) >> 24u;
return res;
}
void main()
{
uint fragInd = gl_GlobalInvocationID.x;
uvec3 fragPos = ExtractLocFromFragment(fragInd);
uint ind = 0u;
int level = int(octreeLevel - 1u);
for (int scale = (1 << (level - 1)), i = 0; i < level; ++i, scale >>= 1)
{
ivec3 idx = ivec3(step(ivec3(scale), fragPos));
ind = (octrees[ind].x & 0x3fffffffu) + (idx.x | (idx.y << 1) | (idx.z << 2));
fragPos &= ~scale;
}
uint brickInd = octrees[ind].y & 0x3fffffffu;
uint blockInd = GetBrickBlock(brickInd);
ivec3 texcoord = ivec3(ExtractBrickLocFromFragment(fragInd));
ivec3 globalCoord = GetBrickCoord(brickInd, texcoord);
octrees[ind].y |= 0x80000000u;
// accumulate color value
uvec3 dstColor = ExtractColorFromFragment(fragInd);
uvec4 srcColor = uvec4(0u);
uint srcVal, dstVal = 0x00000000u;
uint curVal = 0x00000000u;
do
{
srcVal = curVal;
srcColor = DecodeColorFromUint(srcVal);
srcColor = AccumulateColor(srcColor, dstColor);
dstVal = EncodeColorToUint(srcColor);
curVal = imageAtomicCompSwap(brickPool[blockInd].albedos, globalCoord, srcVal, dstVal);
} while (curVal != srcVal);
} 核心代码依旧不复杂，仅仅是在遍历的最后，加上一段写入纹理的操作，使用一个原子操作（imageAtomicCompSwap）来完成。如果追求更高的性能，可以不需要这个操作，直接在对应的纹理坐标上写颜色值就行，虽然会导致flickering。我在插入光线的时候就是直接写纹理数据的（需要达到实时效果，必须减少原子操作）。这样，就算是完成了第一层（LOD0）的数据写入了。接下来，就需要我们手动实现openGL中的glGenerateTextureMipmap函数了。
Isotropic Mipmapping（各项同性过滤）
Interactive Indirect Illumination Using Voxel Cone Tracing论文中，推荐的mipmapping方式，是使用一个3x3的高斯滤波器，对子节点的纹理进行滤波，其中，stride = 2，padding = 1。对父节点中的每个坐标，对子节点的27个值进行采样，采样的结果做加权平均，即可得到父节点的值。注意，高斯滤波器会使用到Brick中的边界值，边界值的传递，我们会在接下来的章节讲述，这里只需要知道Brick中会保存边界值，相邻Brick的边界值相等，就可以了。
const float gaussianWeight[4] = {0.25, 0.125, 0.0625, 0.03125};
const ivec3 brickLocList[27] =
{
ivec3(0, 0, 0),
ivec3(1, 0, 0),
ivec3(2, 0, 0),
ivec3(0, 1, 0),
ivec3(1, 1, 0),
ivec3(2, 1, 0),
ivec3(0, 2, 0),
ivec3(1, 2, 0),
ivec3(2, 2, 0),
ivec3(0, 0, 1),
ivec3(1, 0, 1),
ivec3(2, 0, 1),
ivec3(0, 1, 1),
ivec3(1, 1, 1),
ivec3(2, 1, 1),
ivec3(0, 2, 1),
ivec3(1, 2, 1),
ivec3(2, 2, 1),
ivec3(0, 0, 2),
ivec3(1, 0, 2),
ivec3(2, 0, 2),
ivec3(0, 1, 2),
ivec3(1, 1, 2),
ivec3(2, 1, 2),
ivec3(0, 2, 2),
ivec3(1, 2, 2),
ivec3(2, 2, 2)
};
const float invNormalizedValue = 1.0 / 256.0;
void MipmapIsotropicGaussianFilter(uint nodeInd)
{
for (int ind = 0; ind < 27; ++ind)
{
ivec3 brickLoc = brickLocList[ind];
ivec3 startLoc = 2 * brickLoc - ivec3(1);
vec4 color = vec4(0.0);
float weight = 0.0;
for (int i = 0; i < 27; ++i)
{
ivec3 curBrickLoc = brickLocList[i];
ivec3 bigLoc = startLoc + curBrickLoc;
// 坐标出现溢出
if (bigLoc.x < 0 || bigLoc.y < 0 || bigLoc.z < 0 ||
bigLoc.x > 4 || bigLoc.y > 4 || bigLoc.z > 4)
continue;
int manhattanDist = abs(curBrickLoc.x - 1) + abs(curBrickLoc.y - 1) + abs(curBrickLoc.z - 1);
weight += gaussianWeight[manhattanDist];
ivec3 childLoc = ivec3(step(ivec3(2), bigLoc));
uint childInd = (octrees[nodeInd].x & 0x3fffffffu) + (childLoc.x | (childLoc.y << 1) | (childLoc.z << 2));
// 子节点不存在，跳过
if ((locationList[childInd].y & 0x000000ff) == 0x000000ff)
continue;
ivec3 childBrickLoc = bigLoc - 2 * childLoc;
ivec3 globalCoord = GetBrickCoord((octrees[childInd].y & 0x3fffffffu), childBrickLoc);
vec4 dstColor = imageLoad(brickPool[0].albedos, globalCoord);
color += gaussianWeight[manhattanDist] * dstColor;
}
color /= weight;
color += invNormalizedValue;
ivec3 globalCoord = GetBrickCoord((octrees[nodeInd].y & 0x3fffffffu), brickLoc);
imageStore(brickPool[0].albedos, globalCoord, color);
}
} 实际上，也可以使用另外一种更简单的方法，这个也是openGL实现的Mipmapping算法，就是BoxFilter：
void MipmapIsotropicBoxFilter(uint nodeInd)
{
uint brickInd = octrees[nodeInd].y & 0x3fffffffu;
ivec3 globalCoord = GetBrickCoord(brickInd, ivec3(0));
for (int ind = 0; ind < 8; ++ind)
{
uint childInd = (octrees[nodeInd].x & 0x3fffffffu) + ind;
ivec3 brickLoc = ivec3(ind & 0x1, (ind >> 1) & 0x1, ind >> 2);
if ((locationList[childInd].y & 0x000000ffu) == 0x000000ffu)
imageStore(brickPool[0].albedos, globalCoord + brickLoc, vec4(0.0));
uint childBrickInd = octrees[childInd].y & 0x3fffffffu;
ivec3 targetCoord = GetBrickCoord(childBrickInd, ivec3(0));
vec4 dstColor = vec4(0.0);
for (int i = 0; i < 8; ++i)
{
ivec3 childBrickLoc = ivec3(i & 0x1, (i >> 1) & 0x1, i >> 2);
dstColor += imageLoad(brickPool[0].albedos, targetCoord + childBrickLoc);
}
dstColor *= 0.125;
dstColor += invNormalizedValue;
imageStore(brickPool[0].albedos, globalCoord + brickLoc, dstColor);
}
} 从最终的效果来看，BoxFilter的效果要更好。使用BoxFilter进行滤波的话，就不需要使用Brick中的边界值，只需要对中的值求平均就可以。但是，这两个滤波方法，最终都不会被使用到。因为，SVOGI最终使用的是Anisotropic Mipmapping，也就是各向异性过滤。但是，在自己实现的时候，可以先实现这些简单的滤波来大致观察最终的效果。
Anisotropic Mipmapping（各向异性过滤）
虽然各向同性过滤也能达到不错的效果，实际上github上很多的项目都是使用各向同性过滤来生成纹理的LOD。然而，各项同性过滤会导致一个非常麻烦的问题，就是light leaking。在开头就已经讲过，所谓的Cone Tracing过程，实际上是一个对3D纹理的Ray Marching，每个marching step实际上是对周围场景的一个平均采样，这样会采样大量的空元素，最后采样的alpha值会相对较小。所以这就是light leaking难以解决的原因，各向异性过滤可以有效缓解这个问题。各向异性过滤在原理上也并不复杂。原本，八叉树的每个节点会保存一个3x3x3的纹理块（brick），这个Brick是“无方向”的，也就是mipmapping过程中，在所有方向上的采样频率都是一致。而各向异性过滤下，会为八叉树的每个节点保存6个Brick，每个Brick保存一个方向的结果x，-x，y，-y，z，-z，在mipmapping的过程中，对每个方向上的Brick，针对对应的方向进行alpha blending，其他方向做average blending。听起来好像一头雾水，结合图片就很好理解。比如，一个节点有8个子节点，而每个子节点有一个3x3x3的块，子节点的边界值与邻近节点的对应值是相等的，那么把这8个块连在一块，消除相等的值，就会变成一个5x5x5的块，我们使用一个3x3x3的滤波器，stride=2，padding=1，对这个5x5x5的块进行滤波。
这里使用二维的图像进行展示，实际上3维的情况也差不多。每个父节点像素的mipmapping结果，就是对这个3x3x3的值进行采样的结果。我们之前说过，各向异性过滤就是对对应的方向进行alpha blending，对其他方向做average blending。具体的操作就是对这个3x3x3的Brick，比如这个块是+x方向的，那么此时在x方向上做alpha blending，其他方向上做average blending。我们把这个块再切割成4组值（2维情况下），对这4组值分别在x方向进行alpha blending，在y方向做average blending：
const ivec3 brickLocList[27] =
{
ivec3(0, 0, 0),
ivec3(1, 0, 0),
ivec3(2, 0, 0),
ivec3(0, 1, 0),
ivec3(1, 1, 0),
ivec3(2, 1, 0),
ivec3(0, 2, 0),
ivec3(1, 2, 0),
ivec3(2, 2, 0),
ivec3(0, 0, 1),
ivec3(1, 0, 1),
ivec3(2, 0, 1),
ivec3(0, 1, 1),
ivec3(1, 1, 1),
ivec3(2, 1, 1),
ivec3(0, 2, 1),
ivec3(1, 2, 1),
ivec3(2, 2, 1),
ivec3(0, 0, 2),
ivec3(1, 0, 2),
ivec3(2, 0, 2),
ivec3(0, 1, 2),
ivec3(1, 1, 2),
ivec3(2, 1, 2),
ivec3(0, 2, 2),
ivec3(1, 2, 2),
ivec3(2, 2, 2)
};
const vec3 scaleFactors[27] =
{
vec3(1.0, 1.0, 1.0),
vec3(1.0, 0.5, 0.5),
vec3(1.0, 1.0, 1.0),
vec3(0.5, 1.0, 0.5),
vec3(0.5, 0.5, 0.25),
vec3(0.5, 1.0, 0.5),
vec3(1.0, 1.0, 1.0),
vec3(1.0, 0.5, 0.5),
vec3(1.0, 1.0, 1.0),
vec3(0.5, 0.5, 1.0),
vec3(0.5, 0.25, 0.5),
vec3(0.5, 0.5, 1.0),
vec3(0.25, 0.5, 0.5),
vec3(0.25, 0.25, 0.25),
vec3(0.25, 0.5, 0.5),
vec3(0.5, 0.5, 1.0),
vec3(0.5, 0.25, 0.5),
vec3(0.5, 0.5, 1.0),
vec3(1.0, 1.0, 1.0),
vec3(1.0, 0.5, 0.5),
vec3(1.0, 1.0, 1.0),
vec3(0.5, 1.0, 0.5),
vec3(0.5, 0.5, 0.25),
vec3(0.5, 1.0, 0.5),
vec3(1.0, 1.0, 1.0),
vec3(1.0, 0.5, 0.5),
vec3(1.0, 1.0, 1.0)
};
void MipmapAnisotropicMK2(uint nodeInd, uint idx)
{
uint brickInd = octrees[nodeInd].y & 0x3fffffffu;
ivec3 brickLoc = brickLocList[idx];
ivec3 globalCoord = GetBrickCoord(brickInd, brickLoc);
vec4 colors[8] = {vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0)};
vec4 dstColors[8] = {vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0), vec4(0.0)};
// stride = 2, padding = 1
ivec3 startLoc = 2 * brickLoc - 1;
uint childInd = octrees[nodeInd].x & 0x3fffffffu;
for (int i = 0; i < 8; ++i)
{
ivec3 curLoc = startLoc + ivec3(i & 0x1, (i >> 1) & 0x1, i >> 2);
if (any(lessThan(curLoc, ivec3(0))))
continue;
if (any(greaterThanEqual(curLoc, ivec3(4))))
continue;
ivec3 childLoc = ivec3(step(ivec3(2), curLoc));
uint childOffset = childLoc.x | (childLoc.y << 1) | (childLoc.z << 2);
if ((locationList[childInd + childOffset].y & 0x000000ffu) == 0x000000ffu)
continue;
ivec3 loc = GetBrickCoord(octrees[childInd + childOffset].y & 0x3fffffffu, curLoc - 2 * childLoc);
FetchTexelBlockMK2(0, loc, dstColors);
colors[i] = (dstColors[0] + (1.0 - dstColors[0].a) * dstColors[1] +
dstColors[2] + (1.0 - dstColors[2].a) * dstColors[3] +
dstColors[4] + (1.0 - dstColors[4].a) * dstColors[5] +
dstColors[6] + (1.0 - dstColors[6].a) * dstColors[7]) * 0.25;
}
vec4 pX = (colors[0] + (1.0 - colors[0].a) * colors[1] +
colors[2] + (1.0 - colors[2].a) * colors[3] +
colors[4] + (1.0 - colors[4].a) * colors[5] +
colors[6] + (1.0 - colors[6].a) * colors[7]) * scaleFactors[idx].x;
// repeating another five direction
ivec3 depthOffset = ivec3(0, 0, brickPool[0].depth);
imageStore(brickPool[0].lights, globalCoord, pX);
imageStore(brickPool[0].lights, globalCoord + 1 * depthOffset, nX);
imageStore(brickPool[0].lights, globalCoord + 2 * depthOffset, pY);
imageStore(brickPool[0].lights, globalCoord + 3 * depthOffset, nY);
imageStore(brickPool[0].lights, globalCoord + 4 * depthOffset, pZ);
imageStore(brickPool[0].lights, globalCoord + 5 * depthOffset, nZ);
} 实际上剩下的几个方向有大量的重复代码，这里只贴一个方向（+x）。我使用两个数组保存计算的中间结果，在这两次计算中，我对最后的结果，做了一些调整。在第二次的混合过程中，我把最终结果乘以一个scaleFactors的值，而不是乘以0.25。由于padding = 1，滤波器采样的过程时，在5x5x5以为的纹理，采样的结果是0.0，实际上这里应该采样邻近节点中的对应值，但是为了加速计算，没有访问邻近节点，所以这里用0.0来替代而已。而此时，在第二次计算的时候，如果还是对结果乘以0.25的话，会导致mipmapping的结果过小。实际上，一次的mipmapping的过程是包括节点的mipmapping，以及一次边界值传递的过程来完成。对于边界的值，mipmapping的过程只会保存一个中间的结果，剩下的计算过程会交给Octree Value Transfering来实现。此时，需要乘上一个值，来消减无用值（0.0）的影响，防止边界值过小。比如我们对(0, 0, 0)节点进行mipmapping的时候，其实只会采样到8个节点的值，这时候只需要做一轮的计算，就可以得到最终的结果，而第二轮的计算会退化为对结果乘以0.25，这时候，我们就需要把0.25改为1.0，来获得正确的结果。这个值，在同一个像素的不同方向，也是不同的。比如，对(0, 1, 0)的计算，在+x方向和+y方向是不同的，在+x方向上计算，在第二轮计算的时候，会有两个值剩下，所以scaleFactor=0.5，而在+y方向的计算，在第二轮只会有一个值剩下，所以scaleFactor=1.0。所以每个纹理坐标，需要保存一个3维的向量。听起来似乎很复杂，实际上只要自己做一下这个mipmapping操作，就很容易理解了，对于非中心节点，此时生成的是mipmapping的部分和。
Octree Value Transfering（边界值传递）
在上一章提到，整个八叉树的mipmapping过程，实际上包含两个步骤，一个是对子节点的mipmapping，另一个是边界值的传递。这两个步骤本质上是用一个3x3x3/2x2x2的卷积核，对低一层的像素进行卷积操作的结果，卷积的步长为2（stride=2）。边界值传递（Octree Value Transfering）过程承担着两个重要的任务，第一是把邻近节点的边界值复制到节点的边界，另一个就是对部分mipmapping的像素，完成最后的mipmapping过程。我们先从各项同性过滤（isotropic mipmapping）的角度来看，再从各向异性过滤的角度（anisotropic mipmapping）来观察整个实现。针对不同的mipmapping实现，都有对应不停的transfering算法。我们会一一介绍三个算法，两个是isotropic mipmapping（box filter，gaussian filter），最后一个是anisotropic mipping。
Isotropic Value Transfering
box filter的算法非常简单，只是简单把邻近节点的值复制到该节点的边界像素上，因为全部的mipmapping过程已经在上一步完成了（kernel size=2，stride=2）。这个操作只是为了解决最终vct阶段的时候，纹理采样的问题。
void TransferValueBoxFilter(uint axis, uint nodeInd, ivec3 loc)
{
ivec3 neighborLoc = loc;
uint axis_0 = axis;
neighborLoc += neiTorward[axis_0];
ivec3 neighborSize = textureSize(neighbor, int(hieOffset));
if (any(greaterThanEqual(neighborLoc, neighborSize)))
return;
uint neiNodeInd = texelFetch(neighbor, neighborLoc, int(hieOffset)).r;
if (neiNodeInd == 0xffffffffu)
return;
uint brickInd = octrees[nodeInd].y & 0x3fffffffu;
uint neiBrickInd = octrees[neiNodeInd].y & 0x3fffffffu;
uint blockInd = GetBrickBlock(brickInd);
uint neiBlockInd = GetBrickBlock(neiBrickInd);
uint axis_1 = (axis + 1) % 3;
uint axis_2 = (axis + 2) % 3;
uvec3 remapAxis = uvec3(axis_0, axis_1, axis_2);
ivec3 coord = GetBrickCoord(brickInd, ivec3(0));
ivec3 neiCoord = GetBrickCoord(neiBrickInd, ivec3(0));
ivec3 curBrickLocList[9];
ivec3 neiBrickLocList[9];
for (int ind = 0; ind < 9; ++ind)
{
curBrickLocList[ind] = coord + RemapVec3(ivec3(2, traverseLoc[ind]), remapAxis);
neiBrickLocList[ind] = neiCoord + RemapVec3(ivec3(0, traverseLoc[ind]), remapAxis);
}
vec4 colors[9];
ivec3 depthOffset = ivec3(0, 0, brickPool[0].depth);
FetchTexelBlock(0, neiBrickLocList, colors);
for (int i = 0; i < 9; ++i)
{
imageStore(brickPool[0].albedos, curBrickLocList[i], colors[i]);
}
} 可以看到，只是简单复制邻近节点的值而已。剩下的一个就是gaussian filter，这个tranfering操作也是非常简单，只是在上一个的基础上，对边界的两个值进行一次求平均值而已，这里只贴上不同的部分
vec4 colors[9];
vec4 dstColors[9];
ivec3 depthOffset = ivec3(0, 0, brickPool[0].depth);
FetchTexelBlock(0, neiBrickLocList, colors);
FetchTexelBlock(0, neiBrickLocList, dstColors);
for (int i = 0; i < 9; ++i)
{
vec4 color = 0.5 * colors[i] + 0.5 * dstColors[i];
imageStore(brickPool[0].albedos, curBrickLocList[i], color);
imageStore(brickPool[0].albedos, neiBrickLocList[i], color);
} Anisotropic Transfering
各向异性过滤的做法比之前两个操作都要复杂一些。为此，我们需要再次强调下Anisotropic Mipmapping的核心。Anisotropic Mipmapping需要对六个方向的像素进行mipmapping计算，其中，主要的方向是使用alpha blending，其他方向是使用average blending。作为mipmapping的一部分，值传递步骤也需要进行调整，也就是，主轴方向使用alpha blending，其他轴使用average blending。最好保证，先进行alpha blending，再进行average blending。尽量保证计算结果的正确，这里使用一个数组axisInds，对输入的轴进行重新映射，保证第一轮的传递是alpha blending。
const int axisInds[9] =
{
0, 1, 2,
1, 2, 0,
2, 0, 1
};
void TransferValueAnisotropicMK2(uint axis, uint nodeInd, ivec3 loc, uint dir)
{
ivec3 neighborLoc = loc;
int axis_0 = axisInds[dir * 3 + axis];
neighborLoc += neiTorward[axis_0];
ivec3 neighborSize = textureSize(neighborTex, int(hieOffset));
if (any(greaterThanEqual(neighborLoc, neighborSize)))
return;
uint neiNodeInd = texelFetch(neighborTex, neighborLoc, int(hieOffset)).r;
if (neiNodeInd == 0xffffffffu)
return;
uint brickInd = octrees[nodeInd].y & 0x3fffffffu;
uint neiBrickInd = octrees[neiNodeInd].y & 0x3fffffffu;
uint blockInd = GetBrickBlock(brickInd);
uint neiBlockInd = GetBrickBlock(neiBrickInd);
uint axis_1 = axisInds[dir * 3 + (axis + 1) % 3];
uint axis_2 = axisInds[dir * 3 + (axis + 2) % 3];
uvec3 remapAxis = uvec3(axis_0, axis_1, axis_2);
ivec3 coord = GetBrickCoord(brickInd, ivec3(0));
ivec3 neiCoord = GetBrickCoord(neiBrickInd, ivec3(0));
ivec3 curBrickLocList[9];
ivec3 neiBrickLocList[9];
for (int ind = 0; ind < 9; ++ind)
{
curBrickLocList[ind] = coord + RemapVec3(ivec3(2, traverseLoc[ind]), remapAxis);
neiBrickLocList[ind] = neiCoord + RemapVec3(ivec3(0, traverseLoc[ind]), remapAxis);
}
vec4 srcColors[9];
vec4 colors[9];
ivec3 depthOffset = ivec3(0, 0, brickPool[0].depth);
FetchTexelBlock(dir * 2 + 0, curBrickLocList, srcColors);
FetchTexelBlock(dir * 2 + 0, neiBrickLocList, colors);
vec4 color;
for (int i = 0; i < 9; ++i)
{
// 对主要方向进行alpha blending
if (axis_0 == int(dir))
{
color = srcColors[i] + (1.0 - srcColors[i].a) * colors[i];
imageStore(brickPool[0].lights, curBrickLocList[i] + int(dir * 2 + 0) * depthOffset, color);
imageStore(brickPool[0].lights, neiBrickLocList[i] + int(dir * 2 + 0) * depthOffset, color);
}
// 非主要方向进行average blending
else
{
color = 0.5 * srcColors[i] + 0.5 * colors[i];
imageStore(brickPool[0].lights, curBrickLocList[i] + int(dir * 2 + 0) * depthOffset, color);
imageStore(brickPool[0].lights, neiBrickLocList[i] + int(dir * 2 + 0) * depthOffset, color);
}
}
FetchTexelBlock(dir * 2 + 1, curBrickLocList, srcColors);
FetchTexelBlock(dir * 2 + 1, neiBrickLocList, colors);
for (int i = 0; i < 9; ++i)
{
if (axis_0 == int(dir))
{
color = colors[i] + (1.0 - colors[i].a) * srcColors[i];
imageStore(brickPool[0].lights, curBrickLocList[i] + int(dir * 2 + 1) * depthOffset, color);
imageStore(brickPool[0].lights, neiBrickLocList[i] + int(dir * 2 + 1) * depthOffset, color);
}
else
{
color = 0.5 * srcColors[i] + 0.5 * colors[i];
imageStore(brickPool[0].lights, curBrickLocList[i] + int(dir * 2 + 1) * depthOffset, color);
imageStore(brickPool[0].lights, neiBrickLocList[i] + int(dir * 2 + 1) * depthOffset, color);
}
}
}
// 三个方向的分量，local_size_y = 3
layout (local_size_x = 64, local_size_y = 3, local_size_z = 1) in;
void main()
{
uint threadInd = gl_GlobalInvocationID.x;
uint maxNum = octreeInds[mipmapInd * 2 + 1 + 1];
if (threadInd >= maxNum)
return;
uint begin = octreeInds[mipmapInd * 2 + 1];
uint octreeInd = octreeInds[begin + threadInd];
if ((octrees[octreeInd].y & mipmapMask) == 0u)
return;
ivec3 loc = ExtractLocFromLocationList(octreeInd);
// TransferValueAnisotropic(axis, octreeInd, loc);
TransferValueAnisotropicMK2(0, octreeInd, loc, gl_GlobalInvocationID.y);
memoryBarrier();
TransferValueAnisotropicMK2(1, octreeInd, loc, gl_GlobalInvocationID.y);
memoryBarrier();
TransferValueAnisotropicMK2(2, octreeInd, loc, gl_GlobalInvocationID.y);
} Voxel Cone Tracing
Voxel Cone Tracing本质上是利用3D纹理的线性插值特性，以及八叉树提供的天然mipmapping层级。通过对纹理进行textureLod的过程，实现对一定范围内的值进行采样。svogi的关键其实只是一个八叉树遍历的过程，这个过程在之前的章节中也有介绍，现在，把它使用到tracing中来：
// near：椎体的高
// aperture：椎体半径
vec4 coneTracingAnisotropicMarchingStep(vec3 ori, vec3 dir, float tmin, float tmax, float near, float aperture)
{
float aperNear = (2.0 * aperture) / near;
vec4 color = vec4(0.0);
float width = aperNear * tmin;
vec3 invImgSize = 1.0 / vec3(textureSize(brickPool[0].albedos, 0));
for (; (tmin < tmax) && abs(color.a - 1.0) > EPS; tmin += width)
{
width = aperNear * tmin;
vec3 pos = ori + tmin * dir;
if (IsIllegalPos(pos) == 1)
break;
float lod = max(0.0, log2(width * INV_VOXEL_SIZE));
if (lod > 7.0)
break;
vec3 curPos = 512.0 * pos;
vec3 brickCoord = fract(curPos);
int hiel, hieu;
hiel = int(lod);
hieu = min(hiel + 1, 7);
int scale = 256;
ivec3 octreePos = ivec3(curPos);
ivec3 idx;
uint octreeInd = 0;
vec4 l = vec4(0.0), u = vec4(0.0), a = vec4(0.0);
// 遍历八叉树
for (int level = 8; level > hieu; --level, scale >>= 1)
{
idx = ivec3(step(ivec3(scale), octreePos));
octreeInd = (octrees[octreeInd].x & 0x3fffffffu) + (idx.x | (idx.y << 1) | (idx.z << 2));
if ((locationList[octreeInd].y & 0x000000ffu) == 0x000000ffu)
{
octreeInd = 0xffffffffu;
break;
}
octreePos &= ~scale;
}
// 对节点的值进行采样
if (octreeInd != 0xffffffffu)
{
brickCoord = step(ivec3(scale), octreePos) + brickCoord / float(scale);
TextureBrickAnisotropic(octrees[octreeInd].y & 0x3fffffffu, brickCoord, dir, invImgSize, brickPool[0].lightMap, u);
idx = ivec3(step(ivec3(scale), octreePos));
octreeInd = (octrees[octreeInd].x & 0x3fffffffu) + (idx.x | (idx.y << 1) | (idx.z << 2));
if ((locationList[octreeInd].y & 0x000000ffu) == 0x000000ffu)
{
octreeInd = 0xffffffffu;
break;
}
octreePos &= ~scale;
if (octreeInd != 0xffffffffu)
{
brickCoord = fract(curPos);
brickCoord = step(ivec3(scale), octreePos) + brickCoord / float(scale);
TextureBrickAnisotropic(octrees[octreeInd].y & 0x3fffffffu, brickCoord, dir, invImgSize, brickPool[0].lightMap, l);
}
}
a = mix(l, u, fract(lod));
color += (1.0 - color.a) * a;
}
return color;
} Artifact And Solution
到此为止，已经介绍完了svogi的全部技术细节了。这篇文章不太可能细致地去介绍全部源代码的实现。很多的原理都是来自于后面的几篇reference。这里再讲一下svogi的一个问题（除了light leaking），以及解决方案。svo的一个特点就是，对不存在三角形的节点，是不会生成八叉树节点的。但是，实际上在mipmapping的过程中，存在边界值传递的过程，此时有可能会出现请求邻近节点，但是邻近节点不存在的情况，此时边界值传递失效。不仅如此，在tracing的过程中，如果刚好采样到这个不存在的节点，此时不会采样到任何的值（0.0），然而事实上，该节点是有值的，虽然只会有边界值。此时的tracing结果会出现一个错误，我把它成为tracing撕裂的情况。
可以看到，下方黄色部分，出现了很明显的画面撕裂状况。
解决方法不算复杂，只需要在生成八叉树的时候，对节点周围3个方向的邻近节点都一并生成结果就行了。但是这个做***提高整体的显存开销，在使用的时候，可以具体地看会不会出现这个问题，再采用生成周边节点的方法。
for (uint i = 0; i < 3; ++i)
{
ivec3 tarLoc = coord + DIR[i];
if (tarLoc.x < 0 || tarLoc.y < 0 || tarLoc.z < 0 || tarLoc.x >= scale || tarLoc.y >= scale || tarLoc.z >= scale)
continue;
ivec3 refLoc = tarLoc;
tarLoc <<= (level + 1);
uint neiInd = 0u;
int size = 256;
uint curInd = (octrees[neiInd].x & 0x3fffffffu);
int neiLev = int(octreeLevel - 2);
do
{
size >>= 1u;
int xdiv = tarLoc.x >= size ? 1 : 0;
int ydiv = tarLoc.y >= size ? 1 : 0;
int zdiv = tarLoc.z >= size ? 1 : 0;
neiInd = curInd + uint(xdiv | (ydiv << 1) | (zdiv << 2));
curInd = (octrees[neiInd].x & 0x3fffffffu);
tarLoc.x -= xdiv * size;
tarLoc.y -= ydiv * size;
tarLoc.z -= zdiv * size;
--neiLev;
} while ((neiLev >= 0) && (curInd != 0u));
if (neiLev == -1)
octrees[neiInd].x |= 0x40000000u;
else
octrees[neiInd].x |= 0x80000000u;
imageStore(neighbor, refLoc, uvec4(neiInd));
EncodeLocation(neiInd, refLoc, neiLev + 1);
} Future
todo
reference
- Crassin, Cyril, Fabrice Neyret, Miguel Sainz, Simon Green and Elmar Eisemann: Interactive Indirect Illumination Using Voxel Cone Tracing. In: Computer Graphics Forum, volume 30, pages 1921-1930. Wiley Online Library, 2011.
- NVIDIA Developer. (n.d.). Chapter 42. Conservative Rasterization. [online] Available at: https://developer.nvidia.com/gpugems/gpugems2/part-v-image-oriented-computing/chapter-42-conservative-rasterization.
- Simon (n.d.). Simon’s Tech Blog: Implementing Voxel Cone Tracing. [online] Simon’s Tech Blog. Available at: http://simonstechblog.blogspot.com/2013/01/implementing-voxel-cone-tracing.html [Accessed 20 May 2022].
- Crassin, C. and Green, S. (2012). Octree-Based Sparse Voxelization Using the GPU Hardware Rasterizer. In: OpenGL Insights.