VTK and NumPy – a new take

Ten years ago, I wrote this series of blogs introducing a (then) new way of interfacing VTK and NumPy. In the first blog, I had a teaser that I then expanded on:

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']

Over the last year, we have introduced a simpler and more pythonic interface. As of VTK 9.4, the code above can be expressed as the following.

import vtk
s = (vtk.vtkSphereSource() >> vtk.vtkElevationFilter())()
print(s.point_data.keys())
print(s.point_data['Elevation'])

This code produces the following output.

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

For more details about why line 2 works, see this blog. In summary, we introduced the >> operator for creating pipelines more easily. Furthermore, we introduced the () operator to the pipeline objects (sources, readers, filters etc.) and pipelines for executing and returning the output.

This article is more about lines 3 and 4. Let’s break down line 3. s.point_data is equivalent to s.GetPointData(), just more pythonic. However, if you look at the type of s.point_data, you will notice that it is not vtkPointData, which is what the C++ method returns.

>> print(type(s.point_data))
<class 'vtkmodules.util.data_model.PointData'>
>> print(type(pd).__bases__)
(<class 'vtkmodules.util.data_model.DataSetAttributesBase'>, <class 'vtkmodules.vtkCommonDataModel.vtkPointData'>)

We won’t get into how this is so in this blog. For our purposes, it is sufficient to know that this new class contains new methods such as keys(), values(), and items(). It also implements the [] operator. All the methods used to access the individual arrays can take or return numpy arrays.

elev = s.point_data['Elevation']
print(type(elev))
import numpy
isinstance(elev, numpy.ndarray)

will return

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

Here is a complete example.

import vtk, numpy
s = (vtk.vtkSphereSource() >> vtk.vtkElevationFilter())()
elev = s.point_data['Elevation']
print(numpy.max(elev))
s.point_data['Elev+1'] = elev + 1
print(numpy.max(s.point_data['Elev+1']))

will print

0.5
1.5

I have barely scratched the surface of what is possible in this article. In following articles, I will provide more examples and dig deeper into what is possible.

Credits: Thank you Jaswant Panchumarti, Sebastien Jourdain, and David Gobbi for contributing to this work. This would not be possible without them.