python - 使用 Python 计算表面细胞面积时出现奇怪的 VTK 结果

Question

我有一个关于使用 vtk python 接口(interface)在 VTK 中使用 DeepCopy 和过滤器/算法的问题。我执行以下操作时得到不正确且奇怪的结果:

1. 根据 ConeSource 算法创建圆锥对象 cone
2. 将 cone 深度复制到 cone2
3. 计算圆锥体 cone 和 cone2 的第一个单元格(多边形)的面积。

代码还有一些opt行(opt1、opt2、opt3)。他们改变脚本的输出。

• opt1 暗示导入 vtk_to_numpy，但从未在脚本中使用。
• opt2 改变了cone 的创建方式。通过对新创建的 vtkPolyData 对象应用 DeepCopy 命令(如下面的代码)，或者不进行深度复制，仅引用返回的 vom GetOutput 对象。
• opt3 更改计算和打印第一个单元格面积的次数。

这是最小的工作示例:

import vtk
# from vtk.util.numpy_support import vtk_to_numpy  # opt1

coneSource = vtk.vtkConeSource()
coneSource.Update()

cone = vtk.vtkPolyData()
cone.DeepCopy(coneSource.GetOutput())
# cone = coneSource.GetOutput()  # opt2

cone2 = vtk.vtkPolyData()
cone2.DeepCopy(cone)

N = 80  # opt3
for _ in range(N):
    print cone.GetCell(0).ComputeArea(), cone2.GetCell(0).ComputeArea()

该脚本的结果如下表所示。单元格中的条目是元组，分别表示 cone 和 cone2 的计算区域的输出，其中 ok 是正确结果， (ok) 意味着第一行打印正确，但所有后续计算都不正确。

+--------------------+-----------+-----------+----------+-----------+------------+-----------+
|                    |             DeepCopy             |              GetOutput             |
+--------------------+-----------+-----------+----------+-----------+------------+-----------+
|                    |     1     |    10     |    80    |     1     |     10     |     80    |
+--------------------+-----------+-----------+----------+-----------+------------+-----------+
| w/ numpy_support   |  nan, ok  |  inf, ok  |  ok, ok  |  ok, ok   |  (ok), ok  |  (ok), ok |
| w/o numpy_support  |  inf, ok  |  inf, ok  |  ok, ok  |  ok, inf  |  ok, ok    |  inf, ok  |
+--------------------+-----------+-----------+----------+-----------+------------+-----------+

结果在脚本的后续调用中是可重现的，但在重新启动或较长一段时间后会有所不同。所以你可能会得到不同的结果(尽管我希望你也能体验到一些结果是 nan 或 inf)。

我怀疑我对 VTK 对象做了一些错误的内存操作。如果您能重现我的问题以及我在深复制 vtk 对象或从源创建 vtk 对象期间做错了什么，我对反馈很感兴趣。另外，是否有更好的方法来计算 vtkPolyData 的单元面积，而不是循环遍历所有单元？

我使用VTK-8.0.1、Python2.7并在CentOS7上运行此示例。 VTK和Python都是我自己构建的，所以它们不是存储库中的版本。

Answer 1

我认为 vtkPolygon::ComputeArea() 方法实际上存在错误/意外行为。我看到的问题是 ComputeArea() 的实现如下所示:

double vtkPolygon::ComputeArea()
{
    double normal[3]; //not used, but required for the
                      //following ComputeArea call
    return vtkPolygon::ComputeArea(this->GetPoints(),
                             this->GetNumberOfPoints(),
                             this->GetPointIds()->GetPointer(0),
                             normal);

}

但是，我认为这是“错误的”。 this->GetPoints() 返回仅包含给定多边形的点的数组，因此对于三角形单元来说，这是一个仅包含 3 个点的 0 索引数组。 this->GetPointIds() 是一个带有索引的数组，但如果单元格是大型 polyData 的一部分，则这些不是 this->GetPoints() 的索引数组，而是将其放入包含多数据所有点的全局数组中，而不仅仅是单元格的一个点。即使对于 vtkConeSource 生成的网格的第一个单元格，我发现这是错误的，因为即使第一个单元格使用全局数组中的第一个点(因此访问任何一个都应该没问题)，但它们实际上处于不同的顺序(不知道为什么，但他们是)。

综上所述，如果使用GetPointIds()获取的数组，应该用它来索引“全局”数组。如果使用单元的“本地”vtkPoints，则应从 0 到 N 迭代索引，其中 N 是单元中的点数。 vtkPolygon::ComputeArea() 将这两件事混合在一起，所以我建议此时不要使用它。

如果您将循环内的行替换为使用vtkPolygon::ComputeArea(vtkPoints * p, vtkIdType numPts, vtkIdType * pts, double normal，请验证计算是否正确。 3 ) ) 具有正确的参数？所以像这样的东西(不确定Python语法，抱歉，normal是一个未使用的用于存储法线的向量):

print cone.GetCell(0).ComputeArea(cone.GetPoints(), cone.GetCell(0).GetNumberOfPoints(), cone.GetCell(0).GetPointIds().GetPointer(0), normal), cone2.GetCell(0).ComputeArea(cone2.GetPoints(), cone2.GetCell(0).GetNumberOfPoints(), cone2.GetCell(0).GetPointIds().GetPointer(0), normal)

如果您确实遇到这种情况，我会将问题发布到 vtk bugtracker 上，希望它能在新版本中得到修复。

编辑:您可以使用vtkMeshQuality过滤器来计算面积(以及其他质量指标)。 Here是一个例子，它是用 C++ 编写的，但是针对 python 修改它应该很简单，关键行只是过滤器设置(我对伪 python 的尝试):

qualityFilter = vtk.vtkMeshQuality
qualityFilter.SetInputData(mesh) // mesh in your case is cone/cone2
qualityFilter.SetTriangleQualityMeasureToArea()
qualityFilter.SetQuadQualityMeasureToArea()
qualityFilter.Update()

...以及计算值的检索 - 它们存储在过滤器输出的单元数据数组中:

qualityArray = qualityFilter.GetOutput().GetCellData().GetArray("Quality"))

for i in range(qualityArray->GetNumberOfTuples()):
    print "area of  {0}. element is {1}".format(i, qualityArray->GetValue(i))

就算法复杂度而言，它并不是真的“更好”(如果你希望结果准确，你就无法获得次线性复杂度)，但也许它使用起来更方便。