似乎都需要的头图 细节没怎么扣，有些摩尔纹   前言 每次面试的时候都需要介绍一下自己的工作，有可能是自己的表达能力的关系吧，面试官似乎并不太在意我做的东西，有时候听得云里雾里的。所以自己写一下整个工作的具体实现，既是帮助自己锻炼语言能力，也是求关注，希望得到潜在的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.0vec4 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]);  保守化渲染，三角形的三条边对外扩张后，需要重新计算三角形的三个顶点，保证三角形覆盖到的像素，都能生成片元。在这个阶段，渲染场景需要关闭面裁剪以及深度测试，保证场景的所有三角形都能通过光栅化。保守渲染在自己实现的时候，可能会有一点难度，这里可以使用msaa来代替保守渲染，也能达到一样的效果。  // |-----------|-----------|-------|// 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 = 0for (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  VCTAO的效果  上一章介绍了整个稀疏八叉树的构建过程，这一章，将要介绍的就是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：  在分别计算完这4组值之后，我们此时已经把3x3的值减少为2x2，此时只需要再对这个结果重复一次上面的计算，就可以得到最终的mipmapping结果。对其他方向的操作也是一样的，在所属的放行进行alpha blending，在其他方向做average blending。这里贴出我实现的mipmapping代码，在具体实现上有一些不同：  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 = 3layout (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.