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;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.

The data

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) in a legacy VTK file format. 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.

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

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.

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 and 3 of the already linked examples: These list the main kinds of cell types and their respective indices.

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.

Why not matplotlib?

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.

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. come up with a way to turn all these cell coordinates into surfaces to be plotted with matplotlib, which wouldn't be easy,
2. then realize that 80k surfaces will lead to enormous memory and CPU overhead in matplotlib, finally
3. notice that matplotlib has a 2d renderer, so 3d visualization of complex (read: disjoint, interlocking) surfaces often goes awry.

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

What then?

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.

1. 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 archived wiki page plus these corrections that have changed since (this actually uses a legacy reader, and the code could be modernised quite a bit for more modern input formats, see this newer example):

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()

Note that I've only minimally changed the original wiki version. A more complex example using the legacy reader can be found here. Here's the output after some rotation of the viewport:

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.

2. mayavi.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()

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

Here's how it looks:

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.

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.

3. pyvista

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.

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)

This produces

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.

Answer 2

Let me also draw your attention to

4. vtkplotter

which has a different approach, also built on top of vtk:

from vtkplotter import *

# read the data
ugrid = loadUnStructuredGrid("mesh_mcve.vtk")

# create the outline of the data as polygonal mesh and show it
Mesh(ugrid).c('viridis').alpha(1.0).show()

Plenty of examples can be found here along with API documentation.