VTK和numpy的整合

翻译和修改自：https://blog.kitware.com/improved-vtk-numpy-integration/
英文原文作者：Berk Geveci

Part1

VTK中有一个名为numpy_interface的Python模块，该模块的主要目的是简化VTK和numpy的接口。

import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs# 创建一个默认半径为0.5, 中心位置在(0, 0, 0)的球体
s = vtk.vtkSphereSource()
'''
print(s.GetRadius())
print(s.GetCenter())
'''
# vtkElevationFilter是用于从数据集中生成标量值的过滤器。
e = vtk.vtkElevationFilter()
e.SetInputConnection(s.GetOutputPort())
e.Update()
# 这会将VTK对象包装在一个对象中，该对象可以传递到numpy_interface.algorithms模块中的函数中。
sphere = dsa.WrapDataObject(e.GetOutput())
print(sphere.PointData.keys())
print("Elevation")
print(sphere.PointData['Elevation'])
print("===================================")
print("Normals")
print(sphere.PointData['Normals'])

程序的运行部分结果

<class 'vtk.numpy_interface.dataset_adapter.PolyData'>
['Normals', 'Elevation']
Elevation
[0.5        0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.         0.45048442 0.3117449  0.11126047 0.0.         0.        ]

请注意，我们如何使用其他API访问PointData和Elevation数组。另请注意，当我们打印Elevation数组时，输出看起来与vtkDataArray的输出不同。事实上：

elevation = sphere.PointData['Elevation']
print(type(elevation))
print(isinstance(elevation, numpy.ndarray))

结果如下：

<class 'vtk.numpy_interface.dataset_adapter.VTKArray'>
True

所以vtk数组就是numpy数组？你说的是什么诡计？什么样的魔法能让以下事情成为可能？

sphere.PointData.append(elevation + 1, 'e plus one')
print(algs.max(elevation))
print(algs.max(sphere.PointData['e plus one']))
print(sphere.VTKObject)

部分的输出结果

0.5
1.5
vtkPolyData (0000022A8B1D2E70)....Point Data:....Number Of Arrays: 3Array 0 name = NormalsArray 1 name = ElevationArray 2 name = e plus one

全部都在numpy_interface模块中。它将VTK数据集和数据数组与numpy数组联系在一起，并介绍了许多可以在这些对象上工作的算法。

另外的一个例子：

# 创建用于回归测试的图像
w = vtk.vtkRTAnalyticSource()
# 三角化任何类型的数据集
t = vtk.vtkDataSetTriangleFilter()
t.SetInputConnection(w.GetOutputPort())
t.Update()
ugrid = dsa.WrapDataObject(t.GetOutput())
print(algs.gradient(ugrid.PointData['RTData']))

输出的结果

[[25.46767712  8.78654003  7.28477383][ 6.02292252  8.99845123  7.49668884][ 5.23528767  9.80230141  8.3005352 ]...[-6.43249154 -4.27642226 -8.30053592][-5.19838905 -3.47257614 -7.49668884][13.42047501 -3.26066017 -7.28477287]]

Part2

在这一篇中，将介绍dataset_adapter模块，它是numpy_interface的一部分。该模块旨在简化从Python访问VTK数据集和数组的操作，并提供numpy样式的接口。
使用dataset_adapter模块的第一步是将现有的VTK数据集对象转换为dataset_adapter.VTKObjectWrapper。

s = vtk.vtkSphereSource()
e = vtk.vtkElevationFilter()
e.SetInputConnection(s.GetOutputPort())
e.Update()
sphere = dsa.WrapDataObject(e.GetOutput())
print(sphere)
print(isinstance(sphere, dsa.VTKObjectWrapper))

输出的结果:

<vtk.numpy_interface.dataset_adapter.PolyData object at 0x000002171AB34B70>
True

我们在这里所做的是创建一个dataset_adapter.PolyData类的实例，该实例引用vtkElevationFilter过滤器的输出。我们可以使用VTKObject成员访问基础的VTK对象：

>> print(type(sphere.VTKObject))
<class 'vtkCommonDataModelPython.vtkPolyData'>

WrapDataObject（）函数将为如下的子类返回适当的包装器类:

vtkPolyData
vtkUnstructureGrid
vtkPointSet
vtkDataSet
vtkCompositeDataSet
vtkTable
vtkMolecule
vtkGraph
vtkDataObject

VTKObjectWrapper将VTK方法转发到其VTKObject，以便可以按以下方式直接访问VTK API：

>>> print(sphere.GetNumberOfCells())
96

但是，不能将VTKObjectWrappers作为参数直接传递给VTK方法。

>>> s = vtk.vtkShrinkPolyData()
>>> s.SetInputData(sphere)
Traceback (most recent call last):File "<stdin>", line 1, in <module>
TypeError: SetInputData argument 1: method requires a VTK object
>>> s.SetInputData(sphere.VTKObject)
>>>

数据集属性

到目前为止，很无聊，对吗？我们为VTK数据对象提供了一个包装器，该包装器的行为部分类似于 VTK数据对象。当我们开始寻找如何访问该数据集中包含的字段（数组）时，这会变得更加有趣。

>>> sphere.PointData
<vtk.numpy_interface.dataset_adapter.DataSetAttributes object at 0x000001C36AB7DB38>
>>> sphere.PointData.keys()
['Normals', 'Elevation']
>>> sphere.CellData.keys()
[]
>>> sphere.PointData['Elevation']
VTKArray([0.5       , 0.        , 0.45048442, 0.3117449 , 0.11126047,0.        , 0.        , 0.        , 0.45048442, 0.3117449 ,0.11126047, 0.        , 0.        , 0.        , 0.45048442,0.3117449 , 0.11126047, 0.        , 0.        , 0.        ,0.45048442, 0.3117449 , 0.11126047, 0.        , 0.        ,0.        , 0.45048442, 0.3117449 , 0.11126047, 0.        ,0.        , 0.        , 0.45048442, 0.3117449 , 0.11126047,0.        , 0.        , 0.        , 0.45048442, 0.3117449 ,0.11126047, 0.        , 0.        , 0.        , 0.45048442,0.3117449 , 0.11126047, 0.        , 0.        , 0.        ],dtype=float32)
>>> elevation = sphere.PointData['Elevation']
>>> elevation[:5]
VTKArray([0.5       , 0.        , 0.45048442, 0.3117449 , 0.11126047],dtype=float32)

请注意，这也适用于复合数据集：

>>> mb = vtk.vtkMultiBlockDataSet()
>>> mb.SetNumberOfBlocks(2)
>>> mb.SetBlock(0, sphere.VTKObject)
>>> mb.SetBlock(1, sphere.VTKObject)
>>> mbw = dsa.WrapDataObject(mb)
>>> mbw.PointData
<vtk.numpy_interface.dataset_adapter.CompositeDataSetAttributes object at 0x000001C36ABC4FD0>
>>> mbw.PointData.keys()
['Normals', 'Elevation']
>>> mbw.PointData['Elevation']
<vtk.numpy_interface.dataset_adapter.VTKCompositeDataArray object at 0x000001C36ABD0898>

可以通过这种方式访问PointData，CellData，FieldData，Points（仅vtkPointSet的子类），多边形（仅vtkPolyData）。我们将继续通过此API将访问器添加到更多类型的数组中。

Part3

到目前为止，我已经简要介绍了numpy_interface模块并讨论了数据集接口。最后，我们得到一些更有趣的东西：使用数组，数据集和算法。这是numpy_interface闪耀的地方，并使某些数据分析任务变得更加容易。让我们从一个简单的例子开始。

from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs
import vtkw = vtk.vtkRTAnalyticSource()
w.Update()
image = dsa.WrapDataObject(w.GetOutput())
rtdata = image.PointData['RTData']tets = vtk.vtkDataSetTriangleFilter()
tets.SetInputConnection(w.GetOutputPort())
tets.Update()
ugrid = dsa.WrapDataObject(tets.GetOutput())
rtdata2 = ugrid.PointData['RTData']

在这里，我们创建了两个数据集：图像数据（vtkImageData）和非结构化网格（vtkUnstructuredGrid）。它们本质上表示相同的数据，但非结构化网格是通过对图像数据进行四面体化来创建的。因此，我们希望非结构化网格具有相同的点，但是具有更多的像元（四面体）。

Array API

Numpy_interface数组对象的行为与numpy数组非常相似。实际上，来自vtkDataSet子类的数组是VTKArray的实例，该实例是numpy.ndarray的子类。来自vtkCompositeDataSet和子类的数组不是numpy数组，但是其行为非常相似。

>>> rtdata[0]
60.763466>>> rtdata[-1]
57.113735>>> rtdata[0:10:3]
VTKArray([  60.76346588,   95.53707886,   94.97672272,  108.49817657], dtype=float32)>>> rtdata + 1
VTKArray([ 61.
76346588,  86.87795258,  73.80931091, ...,  68.51051331,44.34006882,  58.1137352 ], dtype=float32)>>> rtdata < 70
VTKArray([ True , False, False, ...,  True,  True,  True], dtype=bool)# 生成矢量场
>>> avector = algs.gradient(rtdata)# 证明avector确实是一个向量
>>> algs.shape(rtdata)
(9261,)>>> algs.shape(avector)
(9261, 3)>>> avector[:, 0]
VTKArray([ 25.69367027,   6.59600449,   5.38400745, ...,  -6.58120966,-5.77147198,  13.19447994])

此示例中需要注意的几件事：

单分量数组始终具有以下形状：（ntuples，）而不是（ntuples，1）
多个组件数组具有以下形状：（ntuples，ncomponents）
张量数组具有以下形状：（ntuples，3，3）
以上内容甚至适用于图像和其他结构化数据。所有数组都具有1维（1个分量数组），2维（多分量数组）或3维（张量数组）。

还有一件很酷的事情。可以使用布尔数组来索引数组。因此，以下代码非常有效：

>>> rtdata[rtdata < 70]
VTKArray([ 60.76346588,  66.75043488,  69.19681549,  50.62128448,64.8801651 ,  57.72655106,  49.75050354,  65.05570221,57.38450241,  69.51113129,  64.24596405,  67.54656982,...,61.18143463,  66.61872864,  55.39360428,  67.51051331,43.34006882,  57.1137352 ], dtype=float32)>>> avector[avector[:,0] > 10]
VTKArray([[ 25.69367027,   9.01253319,   7.51076698],[ 13.1944809 ,   9.01253128,   7.51076508],[ 25.98717642,  -4.49800825,   7.80427408],...,[ 12.9009738 , -16.86548471,  -7.80427504],[ 25.69366837,  -3.48665428,  -7.51076889],[ 13.19447994,  -3.48665524,  -7.51076794]])

算法

只需使用数组API即可完成很多工作。但是，当开始使用numpy_interface.algorithms模块时，事情会变得更加有趣。在前面的示例中对其进行了简要介绍。这里进一步扩展。有关算法的完整列表，请使用help（algs）。以下是一些不言自明的示例：

>>> algs.sin(rtdata)
VTKArray([-0.87873501, -0.86987603, -0.52497   , ..., -0.99943125,-0.59898132,  0.53547275], dtype=float32)>>> algs.min(rtdata)
VTKArray(37.35310363769531)>>> algs.max(avector)
VTKArray(34.781060218811035)>>> algs.max(avector, axis=0)
VTKArray([ 34.78106022,  29.01940918,  18.34743023])>>> algs.max(avector, axis=1)
VTKArray([ 25.69367027,   9.30603981,   9.88350773, ...,  -4.35762835,-3.78016186,  13.19447994])

如果你以前从未使用过axis参数，那么这很简单。当不传aixs值时，该函数将应用于数组的所有值，而无需考虑维数。当axis = 0时，该函数将独立应用于数组的每个组件。当axis = 1时，该功能将独立应用于每个元组。如果不清楚，请尝试一下。以这种方式工作的函数包括sum，min，max，std和var。

另一个有趣且有用的函数是where，它返回发生特定条件的数组的索引。

>>> algs.where(rtdata < 40)
(array([ 420, 9240], dtype=int64),)

对于矢量，如果未定义轴，这还将返回组件索引。

>>> algs.where(avector < -29.7 )
(array([4357, 4797, 4798, 4799, 5239], dtype=int64), array([1, 1, 1, 1, 1], dtype=int64))

到目前为止，我们讨论的所有功能都是由numpy直接提供的。算法模块中包含许多numpy的ufunc。它们都适用于单个数组和复合数据数组。算法还提供了一些功能，其功能与numpy的功能有所不同。这些函数包括叉乘，点乘，矩阵求逆，行列式，特征值，特征向量等。所有这些函数都应用于每个元组，而不是整个数组/矩阵。例如：

>>> algs.determinant(amatrix)
VTKArray([-1221.2732624 ,  -648.48272183,    -3.55133937, ...,28.25770665,  -629.28498775, -1205.81354636])

请注意，以上所有内容仅利用了每个元组的信息，而不依赖于网格。VTK的最大优势之一是其数据模型支持多种网格，其算法通常可在所有这些网格类型上使用。算法模块提供了一些此功能。利用现有的VTK过滤器可以轻松实现其他功能。我在使用梯度之前生成矢量和矩阵。一个例子：

>>> avector = algs.gradient(rtdata)
>>> amatrix = algs.gradient(avector)

这样的函数需要访问包含数组和关联网格的数据集。这是我们在dataset_adapter中使用ndarray子类的原因之一：

>>> rtdata.DataSet
<vtk.numpy_interface.dataset_adapter.DataSet object at 0x000001C36ACA5400>

每个数组都指向包含它的数据集。梯度函数将网格和数组一起使用，numpy也提供了梯度函数。有什么好激动的？好吧，这个：

>>> algs.gradient(rtdata2)
VTKArray([[25.46767712,  8.78654003,  7.28477383],[ 6.02292252,  8.99845123,  7.49668884],[ 5.23528767,  9.80230141,  8.3005352 ],...,[-6.43249154, -4.27642226, -8.30053592],[-5.19838905, -3.47257614, -7.49668884],[13.42047501, -3.26066017, -7.28477287]])
>>> rtdata2.DataSet.GetClassName()
'vtkUnstructuredGrid'

由于VTK的数据模型，需要访问网格的渐变和算法都能工作，无论该网格是均匀网格，曲线网格还是非结构化网格。查看算法模块中的各种功能，以查看使用它可以完成的所有很酷的事情。
所有这些都与复合数据集一起使用，并使用MPI并行进行。我将在以后的博客中介绍有关这些内容的一些具体细节。

Part4

介绍MPI并行分布式编程的，可查看英文原文
https://blog.kitware.com/improved-vtk-numpy-integration-part-4/

Part5

首先定义什么是复合数据集。从类的角度来看，它是vtkCompositeDataSet或其任何子类。从功能的角度来看，这是一种将一组vtkDtaObject（通常是vtkDataSet）集合在一起的方法。最通用的示例是vtkMultiBlockDataSet，它允许创建vtkDataObjects的任意树。另一个示例是vtkOverlappingAMR，它表示Berger-Colella样式的AMR网。这是我们如何创建多块数据集的方法。

>>> import vtk
>>> s = vtk.vtkSphereSource()
>>> s.Update()
>>> c = vtk.vtkConeSource()
>>> c.Update()
>>> mb = vtk.vtkMultiBlockDataSet()
>>> mb.SetBlock(0, s.GetOutput())
>>> mb.SetBlock(1, c.GetOutput())

VTK的许多算法无需更改即可使用复合数据集。例如：

>>> e = vtk.vtkElevationFilter()
>>> e.SetInputData(mb)
>>> e.Update()
>>> mbe = e.GetOutputDataObject(0)
>>> print(mbe.GetClassName())
vtkMultiBlockDataSet

这将输出“ vtkMultiBlockDataSet”。请注意，我使用的是GetOutputDataObject（）而不是GetOutput（）。 GetOutput（）只是一个用SafeDownCast（）包装到算法的预期输出类型的GetOutputDataObject（），在这种情况下为vtkDataSet。因此，即使GetOutputDataObject（）返回实际的复合数据集，GetOutput（）也将返回0。

现在我们有了一个带标量的复合数据集，我们可以使用numpy_interface了。

>>> from vtk.numpy_interface import dataset_adapter as dsa
>>> mbw = dsa.WrapDataObject(mbe)
>>> mbw.PointData.keys()
['Normals', 'Elevation']
>>> elev = mbw.PointData['Elevation']
>>> elev
<vtk.numpy_interface.dataset_adapter.VTKCompositeDataArray object at 0x000001EE12978A58>

请注意，数组类型与我们之前看到的（VTKArray）不同。但是，它仍然以相同的方式工作。

>>> from vtk.numpy_interface import algorithms as algs
>>> algs.max(elev)
0.5
>>> algs.max(elev + 1)
1.5

您可以如下分别访问每个块的数组。

>>> elev.Arrays[0]
VTKArray([0.5       , 0.        , 0.45048442, 0.3117449 , 0.11126047,0.        , 0.        , 0.        , 0.45048442, 0.3117449 ,0.11126047, 0.        , 0.        , 0.        , 0.45048442,0.3117449 , 0.11126047, 0.        , 0.        , 0.        ,0.45048442, 0.3117449 , 0.11126047, 0.        , 0.        ,0.        , 0.45048442, 0.3117449 , 0.11126047, 0.        ,0.        , 0.        , 0.45048442, 0.3117449 , 0.11126047,0.        , 0.        , 0.        , 0.45048442, 0.3117449 ,0.11126047, 0.        , 0.        , 0.        , 0.45048442,0.3117449 , 0.11126047, 0.        , 0.        , 0.        ],dtype=float32)

请注意，索引编制略有不同。

>>> print(elev[0:3])
[VTKArray([0.5       , 0.        , 0.45048442], dtype=float32),
VTKArray([0.       , 0.       , 0.4330127], dtype=float32)]

返回值是一个由2个VTKArrays组成的复合数组。 []运算符仅返回每个数组的前4个值。通常，所有索引操作都适用于复合数组集合中的每个VTKArray。对于where等算法也是如此。

>>> print(algs.where(elev < 0.5))
[(array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49],dtype=int64),), (array([0, 1, 2, 3, 4, 5, 6], dtype=int64),)]

现在，让我们看一下另一个名为Normals的数组。

>>> normals = mbw.PointData['Normals']
>>> normals.Arrays[0]
>>> normals.Arrays[0]
VTKArray([[ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00],[ 0.0000000e+00,  0.0000000e+00, -1.0000000e+00],[ 4.3388373e-01,  0.0000000e+00,  9.0096885e-01],[ 7.8183150e-01,  0.0000000e+00,  6.2348980e-01],[ 9.7492790e-01,  0.0000000e+00,  2.2252093e-01],[ 9.7492790e-01,  0.0000000e+00, -2.2252093e-01],......[ 3.0680212e-01, -3.0680212e-01,  9.0096885e-01],[ 5.5283833e-01, -5.5283833e-01,  6.2348980e-01],[ 6.8937814e-01, -6.8937814e-01,  2.2252093e-01],[ 6.8937814e-01, -6.8937814e-01, -2.2252093e-01],[ 5.5283833e-01, -5.5283833e-01, -6.2348980e-01],[ 3.0680212e-01, -3.0680212e-01, -9.0096885e-01]], dtype=float32)
>>> normals.Arrays[1]
<vtk.numpy_interface.dataset_adapter.VTKNoneArray object at 0x000001CCDE1C0400>

注意第二个数组是VTKNoneArray。这是因为vtkConeSource不产生法线。在不存在数组的地方，我们使用VTKNoneArray作为占位符。这使我们能够在复合数据集的数据集和VTKCompositeDataArray中的数组之间保持一对一的映射。

当许多算法独立地应用于集合中的每个数组时，某些算法是全局的。例如，如上所示，最小和最大。有时获得每个块的答案很有用。为此，您可以使用_per_block算法。

>>> print(algs.max_per_block(elev))
[0.5, 0.4330127]

这些与其他操作一起很好地工作。例如，这是我们如何规范每个块中的高程值。

>>> _max = algs.max_per_block(elev)
>>> _norm = (elev - _min) / (_max - _min)
>>> print(algs.min(_norm))
0.0
>>> print(algs.max(_norm))
1.0
>>>

一旦掌握了这些功能，就应该能够像以前的博客中所述，非常类似于单个阵列使用复合阵列。
关于复合数据集的最后说明。numpy_interface.dataset_adapter提供的复合数据包装器提供了一些方便的功能来遍历复合数据集。这是一个简单的示例：

>>> for ds in mbw:
...   print(type(ds))
...
<class 'vtk.numpy_interface.dataset_adapter.PolyData'>
<class 'vtk.numpy_interface.dataset_adapter.PolyData'>