Improved VTK – numpy integration

July 27, 2014

Recently, we introduced a new Python module called numpy_interface to VTK. The main objective of this module is to make it easier to interface VTK and numpy. This article is the first in a series that introduces this module. Let’s start with a teaser.

import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs
s = vtk.vtkSphereSource()
e = vtk.vtkElevationFilter()
e.SetInputConnection(s.GetOutputPort())
e.Update()
sphere = dsa.WrapDataObject(e.GetOutput())
print sphere.PointData.keys()
print sphere.PointData['Elevation']

This example prints out the following (assuming that you have a relatively new checkout of VTK master from git).

['Normals', '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.        ]

The last three lines are what is new. Note how we used a different API to access the PointData and the Elevation array on the last two lines. Also note that when we printed the Elevation array, the output didn’t look like one from a vtkDataArray. In fact:

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

prints the following.

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

So a VTK array is a numpy array? What kind of trickery is this you say? What kind of magic makes the following possible?

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

Output:

0.5  
1.5  
vtkPolyData (0x7fa20d011c60)  
  ...  
  Point Data:  
    ...  
    Number Of Arrays: 3  
    Array 0 name = Normals  
    Array 1 name = Elevation  
    Array 2 name = e plus one

It is all in the numpy_interface module. It ties VTK datasets and data arrays to numpy arrays and introduces a number of algorithms that can work on these objects. There is quite a bit to this module and I will introduce it piece by piece in upcoming blogs.

Let’s wrap up this blog with one final teaser:

w = vtk.vtkRTAnalyticSource()
t = vtk.vtkDataSetTriangleFilter()
t.SetInputConnection(w.GetOutputPort())
t.Update()
ugrid = dsa.WrapDataObject(t.GetOutput())
print algs.gradient(ugrid.PointData['RTData'])

Output:

[[ 25.46767712   8.78654003   7.28477383]  
 [  6.02292252   8.99845123   7.49668884]  
 [  5.23528767   9.80230141   8.3005352 ]  
 ...,  
 [ -6.43249083  -4.27642155  -8.30053616]  
 [ -5.19838905  -3.47257614  -7.49668884]  
 [ 13.42047501  -3.26066017  -7.28477287]]

Please note that this example is not very easily replicated by using pure numpy. The gradient function returns the gradient of an unstructured grid – a concept that does not exist in numpy. However, the ease-of-use of numpy is there.

Continue on to Part 2.

All posts in this series: Part 1, Part 2, Part 3, Part 4, and Part 5.

3 comments to Improved VTK – numpy integration

  1. I am really glad to see a set of posts covering numpy integration. Very useful information.

    The code blocks in this post are displaying as just an html dump. I see that you mostly fixed the issue in Part 2, but it would be nice if the code showed up clearly here too. Mainly because it is the first post of this 5 part series.

    In Part 2 and I believe other posts, there are still some issues. For example the > character is displaying as a set of four characters. If I type them here, it won’t let me submit the comment.

  2. Hi Berk,

    Thanks for this blog on vtk numpy integration. In most parts it appears very helpful. However, for polydata cell pointIds, the form of the array returned isn’t that useful. For example:

    from vtk.numpy_interface import dataset_adapter as dsa
    polydata_wrapped = dsa.WrapDataObject(polydata)
    polydata_wrapped.GetPolygons()

    which gives:

    VTKArray([ 4, 4161, 0, 4158, 12449, 4, 4158, 1, 4159,
    12449, …, 3570, 12436, 16582], dtype=int64)

    The first number in the array is 4, saying that the cell is a quad and the next 4 numbers [4161, 0, 4158, 12449] are the corresponding points that make up the element. This repeats, with another quad followed by the points [4158, 1, 4159,12449].

    Can you tell me if there is an easy way to turn this into a numpy array or list with the form:
    [[4161, 0, 4158, 12449],
    [4158, 1, 4159,12449],
    …]

    If not, then I think it would be easier not to not use dataset_adapter, but just iterate over the cells and append the cell point to a list.

    Cheers,
    Michael

    1. Hi Michael,

      It is not always easy to work with unstructured data structures like this. If it is all quads, you can achieve what you want with something like this:

      a = numpy.array([4, 1, 2, 3, 4, 4, 5, 6, 7, 8])
      b = a.reshape([a.size // 5, 5])
      b[:, 1:5].flatten()

      There are more complicated ways of dealing with this sort of stuff by converting it to an unstructured grid (using the vtkUnstructuredGrid::SetCells() method), which internally computes an offset array called CellLocations. You can use this array to actually index into the cell array if you have polys that are of different types.

Leave a Reply