如何基于三角形数量动态控制OpenMP并行化三角网格体素化

问题背景 #

需要实现三角网格体素化函数，核心需求：

• 根据三角形数量动态选择是否启用OpenMP并行：数量多时用并行加速，数量少时避免OpenMP启动开销
• 避免if-else导致的循环逻辑代码重复
• 解决非并行场景下#pragma omp critical带来的性能损耗，以及std::unordered_map的线程安全问题

原实现中，if-else重复编写循环逻辑，且匿名函数+临界区的方案会拖慢单线程性能。

可行实现方案 #

1. 用模板消除并行/串行的代码重复 #

利用C++17的if constexpr特性，通过模板参数控制是否启用OpenMP，编译时消除冗余分支，完全避免代码重复：

template<bool UseParallel>
void processTriangles(int NumTriangles, const std::function<void(int)>& processFunc) {
    if constexpr(UseParallel) {
        #pragma omp parallel for num_threads(8)
        for (int i = 0; i < NumTriangles; i++) {
            processFunc(i);
        }
    } else {
        for (int i = 0; i < NumTriangles; i++) {
            processFunc(i);
        }
    }
}

调用方式：根据三角形数量动态选择模板实例，无运行时额外开销

int NumTriangles = vTriangleMesh.getNumTriangles();
if (NumTriangles > 10000) {
    processTriangles<true>(NumTriangles, [&](int i) { /* 三角形处理逻辑 */ });
} else {
    processTriangles<false>(NumTriangles, [&](int i) { /* 三角形处理逻辑 */ });
}

2. 优化临界区与线程安全：线程局部存储+结果合并 #

原代码中#pragma omp critical会导致并行场景下严重的性能瓶颈，且std::unordered_map本身非线程安全。改用线程局部存储+批量合并的方案彻底解决问题：

完整优化后的体素化函数

std::shared_ptr<hive3DModel::CSparseVoxelGrid> voxelizeInSurface(const hive3DModel::CTriangleMesh& vTriangleMesh, double vVoxelSize) {
    _ASSERTE(vVoxelSize > 0);

    Eigen::AlignedBox3d BoundBox;
    for (auto it = vTriangleMesh.vertexBegin(); it != vTriangleMesh.vertexEnd(); it++) {
        BoundBox.extend(it->getPosition());
    }
    Eigen::Vector3d MinBound = BoundBox.min().array() - vVoxelSize / 2;

    auto pVoxelGrid = std::make_shared<hive3DModel::CSparseVoxelGrid>(vVoxelSize, MinBound);
    const Eigen::Vector3d HalfVoxel(vVoxelSize / 2, vVoxelSize / 2, vVoxelSize / 2);
    int NumTriangles = vTriangleMesh.getNumTriangles();
    double InverseVoxelSize = 1 / vVoxelSize;

    // 定义单三角形处理逻辑，返回线程局部体素map
    auto processSingleTriangle = [&](int i) -> std::unordered_map<Eigen::Vector3i, CVoxel, hiveUtils::SEigenHash<Eigen::Vector3i>> {
        std::unordered_map<Eigen::Vector3i, CVoxel, hiveUtils::SEigenHash<Eigen::Vector3i>> localVoxels;
        
        const Eigen::Vector3i& Indexs = vTriangleMesh.getTriangle(i);
        const Eigen::Vector3d& A = vTriangleMesh.getVertex(Indexs[0]).getPosition(),
            B = vTriangleMesh.getVertex(Indexs[1]).getPosition(),
            C = vTriangleMesh.getVertex(Indexs[2]).getPosition();

        Eigen::Vector3d TriangleMinBound = A.cwiseMin(B).cwiseMin(C);
        Eigen::Vector3d TriangleMaxBound = A.cwiseMax(B).cwiseMax(C);

        Eigen::Vector3i MinIndex = Eigen::floor(((TriangleMinBound - MinBound) * InverseVoxelSize).array()).cast<int>();
        Eigen::Vector3i MaxIndex = Eigen::ceil(((TriangleMaxBound - MinBound) * InverseVoxelSize).array()).cast<int>();

        for (int x = MinIndex[0]; x <= MaxIndex[0]; x++) {
            for (int y = MinIndex[1]; y <= MaxIndex[1]; y++) {
                for (int z = MinIndex[2]; z <= MaxIndex[2]; z++) {
                    Eigen::Vector3i Index(x, y, z);
                    Eigen::Vector3d Center = (Index.cast<double>() * vVoxelSize + MinBound) + HalfVoxel;
                    if (hiveUtils::checkTriangleIntersectAABB(A, B, C, Center, HalfVoxel[0])) {
                        localVoxels.insert({Index, CVoxel()});
                    }
                }
            }
        }
        return localVoxels;
    };

    // 动态选择并行/串行处理
    std::vector<std::unordered_map<Eigen::Vector3i, CVoxel, hiveUtils::SEigenHash<Eigen::Vector3i>>> threadLocalVoxels;
    if (NumTriangles > 10000) {
        threadLocalVoxels.resize(8);
        #pragma omp parallel for num_threads(8)
        for (int i = 0; i < NumTriangles; i++) {
            auto local = processSingleTriangle(i);
            auto& threadMap = threadLocalVoxels[omp_get_thread_num()];
            threadMap.insert(local.begin(), local.end());
        }
    } else {
        threadLocalVoxels.emplace_back();
        for (int i = 0; i < NumTriangles; i++) {
            auto local = processSingleTriangle(i);
            threadLocalVoxels[0].insert(local.begin(), local.end());
        }
    }

    // 合并所有线程的局部map到全局体素网格
    for (auto& localMap : threadLocalVoxels) {
        pVoxelGrid->m_Voxels.insert(localMap.begin(), localMap.end());
    }

    return pVoxelGrid;
}

3. 辅助优化：批量插入方法 #

修改CSparseVoxelGrid的addVoxels方法，支持批量插入，减少合并时的函数调用开销：

void addVoxels(const std::unordered_map<Eigen::Vector3i, CVoxel, hiveUtils::SEigenHash<Eigen::Vector3i>>& vVoxels) {
    m_Voxels.insert(vVoxels.begin(), vVoxels.end());
}

内容的提问来源于stack exchange，提问作者user23620257