为什么OpenMP没有并行化vtk IntersectWithLine代码

Posted 2021-05-03

我无法理解为什么如果在循环内部执行Vtk代码，使用OpenMP并行化的for循环不会使用所有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相关的代码？如果没有，可能是什么原因？

答案

它到底是怎么失败的？你有没有尝试过例如打印线程号以查看是否生成了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，所以不是最安全的方式......）