使用python读取和绘制VTK文件数据结构 [英] Reading and plotting VTK file data structure with python

问题描述

我有一个带有points和cells的VTK文件(非结构化网格).

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

我可以导入文件，并使用meshio python包将其读入.

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

如果我键入命令mesh.cells，我会看到一个名为'hexahedron'的字典，其中包含一个由内部列表组成的数组，如下所示:

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]])}

我想在matplotlib中绘制它(我知道ParaView是我一直在使用的替代方法，但目前我也想为此使用matplotlib).无论如何，我在将头包裹在结构上时遇到了麻烦.

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.

每个列表中有8个数据点.

There are 8 data points in each list.

如果运行命令mesh.points，则会得到一个x, y, z坐标列表的数组，这很有意义.但是，对于六面体，列表中是否还有x, y, z坐标?如果有x, y, z坐标列表，那将更有意义，因为这将构成多边形.

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.

随附的是VTK文件，以及在ParaView中是什么样子.谢谢！

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

推荐答案

tl; dr::我认为您不应该尝试使用matplotlib这样做，这将很困难并且不起作用很好.我建议使用专用的vtk库，只需 vtk ，即较高级别的 mayavi.mlab 或我最近获得的最爱

tl;dr: I don't think you should try to use matplotlib for this, and it would be difficult and not work very well. I suggest using a dedicated vtk library, either bare vtk, the higher-level mayavi.mlab or my recently acquired favourite, pyvista. I'll elaborate on all this in detail.

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

First, here is a small, self-contained version of your input data (since the data you linked in the question is both too large and likely to become a broken link sooner or later). I've reduced your data to three rectangular cuboids of varying sizes to approximate your figure.

# 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

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

Let's discuss what this file represents. The header specifies that it's an unstructured grid. This means it can contain points arranged in any kind of arbitrary manner. Basically a bag of points. You can find some explanation about the file format here.

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

The first block, POINTS, contains the 16 float rows, each row corresponds to coordinates of a point in 3d, 16 points in total.

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

The second block, CELLS, defines 3 rows, each row corresponding to a cell (smaller unit, in this case volume) defined in terms of the 0-based indices of the points. The first number (8) indicates the number of vertices in the given cell, the following numbers are the point indices for the corresponding vertices. All three cells in the above example data file consist of 8 vertices, since each cuboid we'll want to draw has 8 vertices. The second number on the CELLS line is the total number of numbers in this block, i.e. 3 * (8+1), i.e. 27.

第三个块CELL_TYPES定义每个3单元格的单元格类型.在这种情况下，它们都是12类型，对应于六面体".内容丰富的人物从已经链接的示例的图2中借来了 : 列出了主要的细胞类型及其各自的索引.

The third block, CELL_TYPES, defines the kind of cell for each of the 3 cells. In this case all of them are type 12, which corresponds to "hexahedrons". An informative figure borrowed from Fig. 2 of the already linked examples: This lists the main kind of cell types and their respective indices.

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

The last block, SCALARS, contains a scalar (number) for each cell according to which it will get coloured later. The scalars 1 through 3 will be mapped onto a colormap to give you the red-to-blue transition seen in your figure.

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

I'm not familiar with meshio but I suspect it gives you access the the aforementioned blocks in the VTK file. The mesh.cells attribute you showed suggests that it recognizes that every cell is a "hexahedron", and lists every cell and their respective 8 vertex indices. The mesh.points attribute is probably an array of shape (n,3), in which case mesh.points[cell_inds, :] gives you the (8,3)-shaped coordinates of a given cell defined by its 8-length array cell_inds.

您如何用matplotlib可视化它?首先，您的实际数据非常庞大，其中包含84480个单元格，尽管从远处看它们看起来很像我上面的示例数据.所以你必须

How would you visualize this with matplotlib? Firstly, your actual data is huge, it contains 84480 cells, even though from afar they look quite like my example data above. So you'd have to

1. 想出一种方法将所有这些单元格坐标转换为要用matplotlib绘制的表面，这并不容易，
2. 然后意识到，最后80k个表面将导致matplotlib中的巨大内存和CPU开销
3. 请注意，matplotlib具有2d渲染器，因此可以对复杂(读取:不相交，互锁)曲面进行3d可视化考虑了所有这些问题，我绝对不会尝试为此使用matplotlib.

   All these things considered, I would definitely not try to use matplotlib for this.

   使用ParaView在幕后使用的东西:VTK！您仍然可以通过低级 vtk 模块或高级( er)级 mayavi.mlab 模块.还有一个mayavi相关的tvtk模块，这是一个中间立场(出于这些目的，它仍然是低级VTK，但具有更python友好的API)，但是我将其保留为练习给读者.

   Use what ParaView uses under the hood: VTK! You can still use the machinery programmatically either via the low-level vtk module or the high(er)-level mayavi.mlab module. There's also the mayavi-related tvtk module which is sort of a middle-ground (it's still low-level VTK for these purposes, but with a more python-friendly API), but I'll leave that as an exercise to the reader.

   使用vtk读取和绘制非结构化网格有点复杂(因为裸vtk总是如此，因为您必须自己组装管道)，但是可以使用

   Reading and plotting an unstructured grid with vtk is a bit complicated (as it always is with bare vtk, since you have to assemble the pipeline yourself), but manageable using this ancient wiki page plus these corrections that have changed since:

   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版本进行了最小改动.这是视口旋转后的输出:

   Note that I've only minimally changed the original wiki version. Here's the output after some rotation of the viewport:

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

   The actual colour depends on the default colormap and the scaling of the scalars. The above default of the vtk module seems to use the jet colormap by default, and it normalizes the scalars so that the values are mapped to the complete colour range.

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

   Personally, I find vtk to be a huge pain to use. It involves a lot of searching, and more often than not digging in the labyrinth of submodules and classes defined in the library. This is why I always try to use vtk via the higher-level functionality of mayavi.mlab instead. This module is especially useful when you're not working with VTK files (like when trying to visualize data that is defined in numpy arrays), but it happens to spare us a lot of work in this case as well, while providing additional functionality. Here's the same visualization using 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上.

   Much less work! We pushed the entire VTK parsing monstrosity onto mayavi, along with the mess of mappers and actors and renderers and ...

   这是它的外观:

   以上是最小的，最省力的可视化，但是您当然可以从此处开始进行任何更改，以使其适合您的需求.您可以更改背景，更改颜色图，以怪异的方式处理数据，然后命名.请注意，这里的颜色与vtk情况相比是相反的，因为默认颜色图或标量到颜色图(查找表)的映射是不同的.您从mlab的高级API越走越多，它就会变得更脏(因为您越来越接近引擎盖下的裸VTK)，但是您通常仍可以使用.

   The above is a minimal, least-effort visualization, but from here of course you can start changing whatever you like to make it fit your needs. You can change the background, change the colormap, manipulate the data in weird ways, you name it. Note that the colours here are reversed compared to the vtk case, since either the default colormap or the mapping of scalars to the colormap (the lookup table) is different. The more you stray from the high-level API of mlab the dirtier it gets (since you're getting closer and closer to the bare VTK under the hood), but you can usually still spare a lot of work and obfuscated code with mayavi.

   最后，mayavi的图形窗口支持所有类型的宝石:交互式修改管道和场景，注释(如坐标轴)，切换正交投影，甚至能够以自动生成的方式记录交互式更改的任何内容python脚本.我绝对建议尝试使用mayavi来实现您要执行的操作.如果您知道使用ParaView会做什么，可以很容易地通过其交互式会话记录功能将其移植到mayavi.

   Finally, mayavi's figure window supports all sorts of gems: interactive modification of the pipeline and scene, annotations such as coordinate axes, toggling orthogonal projections, and even being able to record anything you change interactively in auto-generated python scripts. I would definitely suggest trying to implement what you want to do using mayavi. If you know what you'd do with ParaView, it's fairly easy to port that over to mayavi by making use of its interactive session recording feature.

   最近我指的是pyvista，它是建立在vtk之上的令人愉悦的功能强大的库.尽管它的API需要一些习惯，但是在文档中有很多示例和详尽的API参考 .从库开始有一些学习上的弯路，但是您可能会发现使用它的高级界面和按我的意思"的机制可以提高工作效率.我特别赞赏它的开放源代码性质以及公共问题跟踪程序和响应能力强的核心开发人员.

   I've recently been pointed to pyvista which is a delightfully versatile and powerful library built on top of vtk. Although its API needs some getting used to, there are plenty of examples and an exhaustive API reference in the documentation. There's a bit of a learning curve to get started with the library, but you might find that you're way more productive using its high-level interface and "do what I mean" mechanics. I especially appreciate its open-source nature complete with public issue tracker and responsive core developers.

   那么我们如何使用pyvista读取和绘制网格?去吧:

   So how can we read and plot the grid with pyvista? Here goes:

   import pyvista as pv
   
   # read the data
   grid = pv.read('mesh_mcve.vtk')
   
   # plot the data with an automatically created Plotter
   grid.plot(show_scalar_bar=False, show_axes=False)
   

   这产生

   如您所见，颜色是非常不同的:这是因为pyvista使用matplotlib的感知上统一的 viridis颜色图，非常适合数据可视化！如果您坚持使用怪异的颜色，则可以将cmap='jet'传递给对grid.plot的调用.这里有很多要说的是默认的光照和阴影有何不同，但是我建议仔细阅读文档中有关过滤和绘制数据集的所有选项和方法的信息.

   As you can see the colours are very different: this is because pyvista uses matplotlib's perceptually uniform viridis colormap, which is great for data visualization! If you insist on the weirder colours you can pass cmap='jet' to the call to grid.plot. There's a lot to be said how the default lighting and shading is different here, but I suggest perusing the documentation about all the options and ways to filter and plot datasets.

   这篇关于使用python读取和绘制VTK文件数据结构的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！