为什么 OpenMP 不并行化 vtk IntersectWithLine 代码

Question

如果在循环内执行 Vtk 代码，我很难理解为什么与 OpenMP 并行的 for-loop 不会使用所有 n_threads 线程 (=2x #cores)。具体来说，我想将线/射线与网格相交。我关注了this tutorial

1. 从网格构建 OBB 树
2. 将所有需要的线与网格相交

因为我想并行化它，所以我创建了n_threads 树，这样每个线程都可以使用它自己的树实例：

// Pre-allocate the array
int n_threads = omp_get_max_threads();
trees = std::vector<vtkSmartPointer<vtkOBBTree>>((unsigned int) n_threads);

// Build n_threads OBB trees
#pragma omp parallel for num_threads(n_threads)
for (int t = 0; t < n_threads; ++t)
{
    trees[t] = vtkSmartPointer<vtkOBBTree>::New();
    vtkSmartPointer<vtkPolyData> mesh_clone = vtkSmartPointer<vtkPolyData>::New();
    #pragma omp critical (build_mesh_tree)
    {
        mesh_clone->DeepCopy(mesh);
    }
    trees[t]->SetDataSet(mesh_clone);
    trees[t]->BuildLocator();
}

然后我遍历所有点来计算origin和points中每个点的交集

#pragma omp parallel for num_threads(n_threads)
for (unsigned long i = 0; i < n_points; ++i)
{
    vtkSmartPointer<vtkPoints> intersection_points = vtkSmartPointer<vtkPoints>::New();
    vtkSmartPointer<vtkIdList> cell_ids = vtkSmartPointer<vtkIdList>::New();

    int this_thread = omp_get_thread_num();
    int code = trees[this_thread]->IntersectWithLine(
            origin.data(),       // pointer to the raw data
            points.at(i).data(), // buffer of a Eigen matrix
            intersection_points,
            cell_ids);

    // Do something with intersection_points and cell_ids
}

OpenMP 对于简单的 C++ 代码已经证明可以正常工作。然而，当包裹 Vtk 调用时，它并没有达到它的目的。我想这是因为 Vtk 已经提供了parallelization framework (ref. to the guide)。

如果是这样，您能否解释一下，为什么 OpenMP 无法并行运行与 vtk 相关的代码？如果不是，可能是什么原因？

Answer 1

它究竟是如何失败的？你试过吗？打印线程号以查看是否生成了 n_threads？如果您“只是”在intersection_points 和cell_ids 中得到错误结果，那是因为这些数组在每个线程进入IntersectWithLine 函数时都会被重置（您可以查看实现here，第800 行-808)。

为了解决这个问题，想到的最简单的解决方案是让每个线程都有自己的列表，然后在关键部分将它们连接起来，但更快的方法可能是为每个线程预先分配这些列表的数组，然后读取结果，例如：

vtkSmartPointer<vtkPoints> *inter_points = new vtkSmartPointer<vtkPoints> [n_threads];
vtkSmartPointer<vtkIdList> *inter_cells = new vtkSmartPointer<vtkIdList> [n_threads];
for (unsigned long i = 0; i < n_threads; ++i)
{
    inter_points[i] = vtkSmartPointer<vtkPoints>::New();
    inter_cells[i] = vtkSmartPointer<vtkIdList>::New();
}

#pragma omp parallel for num_threads(n_threads)
for (unsigned long i = 0; i < n_points; ++i)
{
    int this_thread = omp_get_thread_num();
    int code = trees[this_thread]->IntersectWithLine(
            origin.data(), // pointer to the raw data
            points.at(i).data(), // buffer of a Eigen matrix
            inter_points[this_thread],
            inter_cells[this_thread);
}

// now concatenate the results to one array if necessarry

（没有编译，可能有语法错误，只是普通的 C，所以不是最安全的方法...）