如何优化网格法线计算？

Question

我正在为实时网格变形制作应用程序，我需要大量计算法线。现在问题是我通过一些分析发现这段代码占用了最大的 cpu 时间，那么我该如何优化呢？

void Mesh::RecalculateNormals()
{
        for (int i = 0; i < indexCount; i += 3)
        {
            const int ia = indices[i];
            const int ib = indices[i + 1];
            const int ic = indices[i + 2];

            const glm::vec3 e1 = glm::vec3(vert[ia].position) - glm::vec3(vert[ib].position);
            const glm::vec3 e2 = glm::vec3(vert[ic].position) - glm::vec3(vert[ib].position);
            const glm::vec3 no = cross(e1, e2);

            vert[ia].normal += glm::vec4(no, 0.0);
            vert[ib].normal += glm::vec4(no, 0.0);
            vert[ic].normal += glm::vec4(no, 0.0);
        }

        for (int i = 0; i < vertexCount; i++)
            vert[i].normal = glm::vec4(glm::normalize(glm::vec3(vert[i].normal)), 0.0f);

}

此外，在调用此函数之前，我必须循环遍历所有顶点并通过将法线设置为 vec3(0) 来清除先前的法线。

如何加快速度？有没有更好的算法？还是 GLM 很慢？

Answer 1

优化此代码的主要方法有 3 种：使用 SIMD 指令、使用多个 线程 和处理 内存访问模式。它们都不是灵丹妙药。

在您的情况下，使用 SIMD 指令并非易事，因为在第一个循环中内存中基于 indices 的间接数据相关读取。话虽如此，最近的 SIMD 指令集（如 x86-64 处理器上的 AVX-2/AVX-512 和 ARM 处理器上的 SVE）提供了执行收集负载的指令。一旦将值加载到 SIMD 寄存器中，您就可以非常快速地计算叉积。问题是第一个循环中最后一个基于 indices 的间接数据相关存储需要 scatter 存储指令，该指令仅在支持 AVX-512 (x86-64) 和 SVE (手臂）。您可以从 SIMD 寄存器中提取值以便串行存储它们，但这肯定会非常低效。希望第二个循环可以更容易地向量化，但您当然需要以对 SIMD 更友好的方式重新排序正常数据结构（请参阅 AoS vs SoA 和 Data-oriented design）。最后，我希望只要不使用分散指令，第一个循环的 SIMD 实现不会快得多，而第二个循环肯定会快得多。实际上，我希望即使使用收集/分散指令，第一个循环的 SIMD 实现也不会快很多，因为这些指令的实现效率往往很低，而且我希望你的代码的第一个循环会导致很多 缓存未命中（请参阅下一节）。

使用多线程也不是小事。事实上，虽然第一个循环的第一部分可以简单地并行化，但执行法线累加的第二部分不能。一种解决方案是使用原子添加。另一种解决方案是使用每个线程的法线向量数组。更快的解决方案是使用分区方法将网格分割成独立的部分（每个线程都可以安全地执行累积，除了可能共享的vect 项目）。请注意，有效地划分网格以平衡许多线程之间的工作远非简单，AFAIK 已成为许多过去研究论文的主题。最佳解决方案很大程度上取决于线程数、内存中整体vert 缓冲区的大小以及您的性能/复杂性要求。第二个循环可以简单地并行化。并行化循环的最简单解决方案是使用 OpenMP（尽管有效的实现可能相当复杂，但很少有 pragma 足以并行化循环）。

关于输入数据，第一个循环可能非常低效。实际上，ia、ib 和ic 可以引用非常不同的索引，从而导致可预测的vert[...] 在内存中加载/存储。如果结构很大，加载/存储将导致缓存未命中。由于 RAM 的巨大延迟，数据结构适合 RAM，因此此类缓存未命中可能会非常缓慢。解决此问题的最佳解决方案是改善内存访问的局部性。重新排序 vert 项目和/或 indices 可以极大地提高局部性和性能。同样，可以使用分区方法来做到这一点。一个天真的开始是对vert 进行排序，以便在更新ib 和ic 索引时对ia 进行排序，以便它们仍然有效。这可以使用键值 arg-sort 来计算。如果网格是动态的，则问题变得非常复杂，难以有效解决。八叉树和 k-D 树有助于提高计算的局部性，而不会引入（太大）开销。

请注意，您可能不需要重新计算有关先前操作的所有法线。如果是这样，您可以跟踪需要重新计算的那个，并且只执行增量更新。

最后，请检查编译器优化是否启用。此外，请注意，您可以使用（不安全的）fash-math 标志来改进代码的自动矢量化。您还应该检查使用的 SIMD 指令集以及编译器是否内联了 glm 调用。

Answer 2

虽然 Jérôme Richard 的建议要好得多，但到目前为止，实施对我来说很复杂的一切。所以我尝试了一些基本的优化，现在代码快了大约 5 到 6 倍！

这是新代码：

#define VEC3_SUB(a, b, out) out.x = a.x - b.x; \
                out.y = a.y - b.y; \
                out.z = a.z - b.z;

#define VEC3_ADD(a, b, out) out.x = a.x + b.x; \
                out.y = a.y + b.y; \
                out.z = a.z + b.z;

float inline __declspec (naked) __fastcall asm_sqrt(float n)
{
    _asm fld dword ptr [esp+4]
    _asm fsqrt
    _asm ret 4
} 

#define VEC3_NORMALIZE(v, out)  \
{               \
                                 \
    float tempLength = ( (v.x) * (v.x) ) + ( (v.y) * (v.y) ) +( (v.z) * (v.z) ); \
    float lengthSqrauedI = 1.0f / asm_sqrt(tempLength); \
    out.x = (v.x) * lengthSqrauedI;              \
    out.y = (v.y) * lengthSqrauedI;              \
    out.z = (v.z) * lengthSqrauedI;              \
}

#define VEC3_CROSS(a, b, out) \
{                             \
    out.x = ( (a.y) * (b.z) ) - ( (a.z) * (b.y) ); \
    out.y = ( (a.z) * (b.x) ) - ( (a.x) * (b.z) ); \
    out.z = ( (a.x) * (b.y) ) - ( (a.y) * (b.x) ); \
}

void Mesh::RecalculateNormals()
{
        glm::vec3 e1;
        glm::vec3 e2;   
        glm::vec3 no;

        int iabc[3];
        for (int i = 0; i < indexCount; i += 3)
        {
            iabc[0] = indices[i];
            iabc[1] = indices[i + 1];
            iabc[2] = indices[i + 2];

            glm::vec4& tmp4a = vert[iabc[0]].position;
            glm::vec4& tmp4b = vert[iabc[1]].position;
            glm::vec4& tmp4c = vert[iabc[2]].position;
                        
            VEC3_SUB(tmp4a, tmp4b, e1);
            VEC3_SUB(tmp4c, tmp4b, e2);

            VEC3_CROSS(e1, e2, no);

            VEC3_ADD(vert[iabc[0]].normal, no, vert[iabc[0]].normal);
            VEC3_ADD(vert[iabc[1]].normal, no, vert[iabc[1]].normal);
            VEC3_ADD(vert[iabc[2]].normal, no, vert[iabc[2]].normal);
        }

        for (int i = 0; i < vertexCount; i++)
        {
            VEC3_NORMALIZE(vert[i].normal, vert[i].normal);
        }

}