Reading VTK files in Python via python-vtk

by Forrest Sheng Bao http://fsbao.net

Surprisingly, I noticed that there isn't a good document covering frequently-used functions of the Python module vtk (provided via python-vtk on Debian/Ubuntu Linux systems). So I decide to write a very simple one here, covering all functions that I have used.

Note:

• This tutorial is for Python 2.X though slight changes can make it work with Python 3.X.
• I assume you are familiar with VTK data format, thus you know what header, DATASET and POINT_DATA are.
• Notations like In [123] or Out[123] are prompts in iPython (an interactive Python shell). They show a line of code and its output/effect, respectively. They should NOT appear in your code. And, when you use other Python interpreter/shell, you may not see it.
• I did not use print function in iPython when showing something. But when you write a Python program, you need print to display.

Step 1: Set up the reader.

import vtk
reader = vtk.vtkDataSetReader()
reader.SetFileName("lh.sulc.fundi.from.pits.pial.vtk")
reader.ReadAllScalarsOn()  # Activate the reading of all scalars
reader.Update()

The reader is the top level object we use to access a VTK file. For example, you can know the header of a VTK file by

In [119]: reader.GetHeader()
Out[119]: 'Created by Mindboggle'

Step 2: Get data in DATASET block.

Again, reader is the top level. To get data in DATASET block, we need to call a method on reader.

data=reader.GetOutput()

We can know how many (not ``much'' here) data of each type are in the DATASET by the function GetNumberOf{VTK_Data_Type}(). For example, to know how may Points are there, we use:

In [11]: data.GetNumberOfPoints()
Out[11]: 128895

I used tab-completion function of iPython to find out all GetNumberOf{VTK_Data_Type}() for data:

data.GetNumberOfCells   data.GetNumberOfPolys
data.GetNumberOfLines   data.GetNumberOfStrips
data.GetNumberOfPieces  data.GetNumberOfVerts
data.GetNumberOfPoints  

After knowing the size of data of each type, we can access them using function Get{VTK_DATA_Type}(). For example, to get the first point's coordinate, we can do this:

In [136]: data.GetPoint(0)
Out[136]: (-11.605026245117188, -97.47259521484375, 3.8222298622131348)

Step 3: get data in POINT_DATA/CELL_DATA block.

Data in POINT_DATA and CELL_DATA can be loaded via function GetPointData() and GetCellData(). For example,

d=data.GetPointData()

We will use loading a scalar array from POINT_DATA as an example below.
Vectors and tensors can be accessed in similar way by different functions.

All scalar arrays in a file can be viewed by:

In [140]: reader.GetNumberOfScalarsInFile() # get number of scalars
Out[140]: 5
In [141]: reader.GetScalarsNameInFile(1) # get scalar name string
Out[140]: curv

A scalar can be accessed by its name:

In [38]: array=d.GetArray('curv') # 'curv' is the scalar name
In [50]: array.GetValue(260957-260761+1)
Out[50]: 0.043181419372558594

In the end, you should see the following variables:

Variable   Type         Data/Info
---------------------------------
curv       vtkobject    vtkFloatArray (0x9e5cd38)<...>)\n  Array: 0x9f1a9a0\n\n
d          vtkobject    vtkPointData (0x9789b50)\<...>  PedigreeIds: (none)\n\n
data       vtkobject    vtkPolyData (0x9e6c650)\n<...>: 0\n  Ghost Level: 0\n\n
reader     vtkobject    vtkDataSetReader (0x9c431<...> InputStringLength: 0\n\n
vtk        module       < module 'vtk' from '/usr/<...>hon2.6/vtk/__init__.pyc'>

The end. Comments and questions are welcomed.

References:

• http://www.vtk.org/pipermail/vtkusers/2010-June/109492.html
• http://stackoverflow.com/questions/6684306/how-can-i-read-a-vtk-file-into-a-python-datastructure