Reading and plotting VTK file data structure with python

Question

I have a VTK file (unstructured grid) with points and cells.

I can import the file and read it into using the meshio python package.

If I type the command mesh.cells I see a dictionary called 'hexahedron' with an array made up of lists inside like this:

{'hexahedron': array([[  0, 162, 185, ..., 163, 186,  23],
        [162, 329, 351, ..., 330, 352, 186],
        [329, 491, 514, ..., 492, 515, 352],
        ...,
        [483, 583, 600, ..., 584, 601, 490],
        [583, 650, 656, ..., 651, 657, 601],
        [650, 746, 762, ..., 747, 763, 657]])}

I would like to plot this in matplotlib (I know ParaView is an alternative, which I've been using, but I would also like to use matplotlib for this at the moment). Anyways, I'm having trouble wrapping my head around the structure.

There are 8 data points in each list.

If I run the command mesh.points I get an array of lists of x, y, z coordinates, which makes sense. However, with the hexahedron, are there also x, y, z coordinates in the list? It would make more sense if there were lists of x, y, z coordinates, as that would make up polygons.

I've seen this thread, but I'm still stuck on understanding this.

Attached is the VTK file, as well as what it looks like in ParaView. Thanks!

Answer 1

tl;博士：我不认为您应该为此使用matplotlib，这将很困难并且不能很好地工作。我建议使用专用的vtk库，无论是裸库vtk还是更高级别的库mayavi.mlab。我将详细说明所有这些。

数据

首先，这是输入数据的一个很小的独立版本（因为您在问题中链接的数据太大，并且很可能迟早会成为断开的链接）。我已将您的数据缩小为三个大小不同的矩形长方体，以近似您的身材。

# vtk DataFile Version 3.1
MCVE VTK file
ASCII
DATASET UNSTRUCTURED_GRID
POINTS      16 float
 0.   0.   0.
 0.   0.   3.
 0.   2.   0.
 0.   2.   3.
 4.   0.   0.
 4.   0.   3.
 4.   2.   0.
 4.   2.   3.
 5.   0.   0.
 5.   0.   3.
 5.   2.   0.
 5.   2.   3.
13.   0.   0.
13.   0.   3.
13.   2.   0.
13.   2.   3.

CELLS        3     27
 8    0   1   3   2   4   5   7   6
 8    4   5   7   6   8   9  11  10
 8    8   9  11  10  12  13  15  14

CELL_TYPES        3
          12          12          12

CELL_DATA        3
SCALARS elem_val float
LOOKUP_TABLE default
 1
 2
 3

让我们讨论这个文件代表什么。标头指定它是非结构化网格。这意味着它可以包含以任何一种任意方式排列的点。基本上是一袋积分。您可以在此处找到有关文件格式的一些说明。

第一块POINTS包含16 float行，每行对应3d中一个点的坐标，总共16个点。

第二个块CELLS定义3行，每行对应于一个基于点的从0开始的索引所定义的单元格（在这种情况下为体积较小的单元）。第一个数字（8）表示给定像元中顶点的数量，随后的数字是对应顶点的点索引。上面的示例数据文件中的所有三个单元都包含8个顶点，因为我们要绘制的每个长方体都有8个顶点。该CELLS行上的第二个数字是此块中的总数3 * (8+1)，即27。

第三块CELL_TYPES定义每个单元的3单元类型。在这种情况下，它们都是type 12，对应于“ hexahedrons”。从已经链接的示例的图2中借来的信息量很大： 这列出了主要的细胞类型及其各自的索引。

最后一个块SCALARS包含每个单元格的标量（数字），以后将根据该标量进行着色。标量1通过3将映射到颜色图上，从而为您提供从图中看到的红色到蓝色的过渡。

为什么不使用matplotlib？

我不熟悉，meshio但我怀疑它使您可以访问VTK文件中的上述块。mesh.cells您显示的属性表明它识别出每个单元格都是“六面体”，并列出了每个单元格及其各自的8个顶点索引。该mesh.points属性可能是一个形状数组(n,3)，在这种情况下，mesh.points[cell_inds, :]将为您(8,3)提供由其8长度数组定义的给定单元格的形坐标cell_inds。

您如何用matplotlib可视化？首先，您的实际数据非常庞大，其中包含84480个单元，即使从远处看它们看起来与我上面的示例数据非常相似。所以你必须

1. 想出一种方法将所有这些单元格坐标转换成要用matplotlib绘制的表面，这并不容易，
2. 然后意识到80k曲面将在matplotlib中导致巨大的内存和CPU开销，最后
3. 请注意，matplotlib具有2d渲染器，因此复杂（读取，不相交，互锁）曲面的3d可视化通常会出错。

考虑到所有这些因素，我绝对不会尝试为此使用matplotlib。

然后怎样呢？

使用ParaView在幕后使用的东西：VTK！您仍然可以通过低级别vtk模块或高级别模块以编程方式使用机械mayavi.mlab。还有一个mayavi-related tvtk模块，它是一个中间立场（出于这些目的，它仍然是低级VTK，但是具有更Python友好的API），但我将其留给读者练习。

1。 vtk

使用vtk读取和绘制非结构化网格有点复杂（因为裸vtk总是如此，因为您必须自己组装管道），但是可以使用此古老的Wiki页面以及自以下以来已更改的更正内容进行管理：

from vtk import (vtkUnstructuredGridReader, vtkDataSetMapper, vtkActor,
                 vtkRenderer, vtkRenderWindow, vtkRenderWindowInteractor)

file_name = "mesh_mcve.vtk"  # minimal example vtk file

# Read the source file.
reader = vtkUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()  # Needed because of GetScalarRange
output = reader.GetOutput()
output_port = reader.GetOutputPort()
scalar_range = output.GetScalarRange()

# Create the mapper that corresponds the objects of the vtk file
# into graphics elements
mapper = vtkDataSetMapper()
mapper.SetInputConnection(output_port)
mapper.SetScalarRange(scalar_range)

# Create the Actor
actor = vtkActor()
actor.SetMapper(mapper)

# Create the Renderer
renderer = vtkRenderer()
renderer.AddActor(actor)
renderer.SetBackground(1, 1, 1) # Set background to white

# Create the RendererWindow
renderer_window = vtkRenderWindow()
renderer_window.AddRenderer(renderer)

# Create the RendererWindowInteractor and display the vtk_file
interactor = vtkRenderWindowInteractor()
interactor.SetRenderWindow(renderer_window)
interactor.Initialize()
interactor.Start()

请注意，我仅对原始Wiki版本进行了最小程度的更改。这是视口旋转后的输出：

实际颜色取决于默认颜色图和标量的缩放比例。上面的默认vtk模块似乎默认情况下使用jet颜色图，并且对标量进行了规范化，以便将值映射到整个颜色范围。

2。 mayavi.mlab

就个人而言，我觉得vtk使用起来非常痛苦。它涉及大量的搜索，并且更多地是在库中定义的子模块和类的迷宫中进行挖掘。这就是为什么我总是尝试vtk通过更高级别的功能来使用它的原因mayavi.mlab。当您不使用VTK文件时（例如，当尝试可视化numpy数组中定义的数据时），该模块特别有用，但是在这种情况下，它在提供其他功能的同时也为我们节省了很多工作。这是使用的相同可视化mlab：

from mayavi import mlab
from mayavi.modules.surface import Surface

file_name = "mesh_mcve.vtk"  # minimal example vtk file

# create a new figure, grab the engine that's created with it
fig = mlab.figure()
engine = mlab.get_engine()

# open the vtk file, let mayavi figure it all out
vtk_file_reader = engine.open(file_name)

# plot surface corresponding to the data
surface = Surface()
engine.add_filter(surface, vtk_file_reader)

# block until figure is closed
mlab.show()

工作少得多！我们将整个VTK解析怪兽推到了mayavi一起，以及一堆乱七八糟的映射器，演员和渲染器以及...

外观如下：

上面是最小的，省力的可视化，但是您当然可以从这里开始进行任何更改，以使其适合您的需求。您可以更改背景，更改颜色图，以奇怪的方式处理数据，然后命名。请注意，与vtk默认情况相比，此处的颜色是相反的，因为默认颜色图或标量到颜色图（查找表）的映射是不同的。从mlab的高级API中流失的次数越多，它就会变得越脏（因为您越来越接近幕后的裸VTK），但是通常您仍然可以节省大量工作和混淆代码mayavi。

最后，mayavi的图形窗口支持各种宝石：交互式修改管道和场景，注释（例如坐标轴），切换正交投影，甚至还可以在自动生成的python脚本中交互式记录您所做的任何更改。我绝对建议尝试使用mayavi实现您想做的事情。如果您知道使用ParaView会做什么，可以很容易地mayavi通过使用其交互式会话记录功能将其移植到。