为什么不的OpenMP并行VTK IntersectWithLine代码

Question

我很难理解为什么for循环，使用OpenMP并行化，不若VTK代码在循环内执行使用所有n_threads线程（= 2X #cores）。具体而言，我想用网格相交线条/光线。我通过

1. 建设OBB树从网格
2. 所有需要的线相交，其中网

因为我想并行它遵循this tutorial，我创建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，所以不是最安全的方式......）