如何优化网格法线计算？

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

优化此代码的三种主要方法：使用 SIMD 指令，使用多个线程并在 内存访问模式。 None 个是 silver-bullet.

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

使用多线程也很重要。事实上，虽然第一个循环的第一部分可以简单地并行化，但执行法线累加的第二部分却不能。一种解决方案是使用原子添加。另一种解决方案是使用 per-thread 法向量数组。更快的解决方案是使用分区方法，将网格分割成独立的部分（每个线程都可以安全地执行累加，possibly-shared vect 项除外）。请注意，有效地划分网格以平衡多个线程之间的工作远非简单，据我所知，它一直是许多过去研究论文的主题。最佳解决方案很大程度上取决于线程数、内存中 vert 缓冲区的总大小以及您的 performance/complexity 要求。第二个循环可以简单地并行化。并行化循环的最简单解决方案是使用 OpenMP（很少 pragma 足以并行化循环，尽管有效的实现可能非常复杂）。

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

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

最后，请检查是否启用了编译器优化。此外，请注意，您可以使用（不安全的）fash-math 标志来改进代码的 auto-vectorization。您还应该检查使用的 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);
        }

}

如何优化网格法线计算？

How to optimize mesh normals calculation?

c++

algorithm

optimization

mesh

normals