Improved VTK - numpy integration (part 5)

Welcome to my last blog in the series where we to discover VTK's numpy_interface module. If you are not familiar with this module, I recommend checking out my previous blogs on it ([1], [2], [3]). In this blog, I will talk about how one can work with composite datasets and arrays using this module.

Let's start with defining what a composite dataset is. From a class point of view, it is vtkCompositeDataSet or any of its subclasses. From a functionality point of view, it is a way of collecting together a set of vtkDataObjects (usually vtkDataSets). The most generic example is vtkMultiBlockDataSet which allows the creation of an arbitrary tree of vtkDataObjects. Another example is vtkOverlappingAMR which represent a Berger-Colella style AMR meshes. Here is how we can create a multi-block dataset.

>>> import vtk
>>> s = vtk.vtkSphereSource()
>>> s.Update()
>>> c = vtk.vtkConeSource()
>>> c.Update()
>>> mb = vtk.vtkMultiBlockDataSet()
>>> mb.SetBlock(0, s.GetOutput())
>>> mb.SetBlock(1, c.GetOutput())

Many of VTK's algorithms work with composite datasets without any change. For example:

>>> e = vtk.vtkElevationFilter()
>>> e.SetInputData(mb)
>>> e.Update()
>>> mbe = e.GetOutputDataObject(0)
>>> print mbe.GetClassName()

This will output 'vtkMultiBlockDataSet'. Note that I used GetOutputDataObject() rather than GetOutput(). GetOutput() is simply a GetOutputDataObject() wrapped with a SafeDownCast() to the expected output type of the algorithm - which in this case is a vtkDataSet. So GetOutput() will return 0 even when GetOutputDataObject() returns an actual composite dataset.

Now that we have a composite dataset with a scalar, we can use numpy_interface.

>>> from vtk.numpy_interface import dataset_adapter as dsa
>>> mbw = dsa.WrapDataObject(mbe)
>>> mbw.PointData.keys()
['Normals', 'Elevation']
>>> elev = mbw.PointData['Elevation']
>>> elev
<vtk.numpy_interface.dataset_adapter.VTKCompositeDataArray at 0x1189ee410>

Note that the array type is different than we have seen previously (VTKArray). However, it still works the same way.

>>> from vtk.numpy_interface import algorithms as algs
>>> algs.max(elev)
0.5
>>> algs.max(elev + 1)
1.5

You can individually access the arrays of each block as follows.

>>> elev.Arrays[0]
VTKArray([ 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.        ], dtype=float32)

Note that indexing is slightly different.

>>> print elev[0:3]
[VTKArray([ 0.5,  0.,  0.45048442], dtype=float32),
 VTKArray([ 0.,  0.,  0.43301269], dtype=float32)]

The return value is a composite array consisting of 2 VTKArrays. The [] operator simply returned the first 4 values of each array. In general, all indexing operations apply to each VTKArray in the composite array collection. Similarly for algorithms such as where.

>>> print algs.where(elev < 0.5)
[(array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
       35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]),),
       (array([0, 1, 2, 3, 4, 5, 6]),)]

Now, let's look at the other array called Normals.

>>> normals = mbw.PointData['Normals']
>>> normals.Arrays[0]
VTKArray([[  0.00000000e+00,   0.00000000e+00,   1.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,  -1.00000000e+00],
       [  4.33883727e-01,   0.00000000e+00,   9.00968850e-01],
       [  7.81831503e-01,   0.00000000e+00,   6.23489797e-01],
       [  9.74927902e-01,   0.00000000e+00,   2.22520933e-01],
       ...
       [  6.89378142e-01,  -6.89378142e-01,   2.22520933e-01],
       [  6.89378142e-01,  -6.89378142e-01,  -2.22520933e-01],
       [  5.52838326e-01,  -5.52838326e-01,  -6.23489797e-01],
       [  3.06802124e-01,  -3.06802124e-01,  -9.00968850e-01]], dtype=float32)
>>> normals.Arrays[1]
<vtk.numpy_interface.dataset_adapter.VTKNoneArray at 0x1189e7790>

Notice how the second arrays is a VTKNoneArray. This is because vtkConeSource does not produce normals. Where an array does not exist, we use a VTKNoneArray as placeholder. This allows us to maintain a one-to-one mapping between datasets of a composite dataset and the arrays in the VTKCompositeDataArray. It also allows us to keep algorithms working in parallel without a lot of specialized code (see my previous blog).

Where many of the algorithms apply independently to each array in a collection, some algorithms are global. For example, min and max as we demonstrated above. It is sometimes useful to get per block answers. For this, you can use _per_block algorithms.

>>> print algs.max_per_block(elev)
[VTKArray(0.5, dtype=float32), VTKArray(0.4330126941204071, dtype=float32)]

These work very nicely together with other operations. For example, here is how we can normalize the elevation values in each block.

>>> _min = algs.min_per_block(elev)
>>> _max = algs.max_per_block(elev)
>>> _norm = (elev - _min) / (_max - _min)
>>> print algs.min(_norm)
0.0
>>> print algs.max(_norm)
1.0

Once you grasp these features, you should be able to use composite array very similarly to single arrays as described in previous blogs.

A final note on composite datasets. The composite data wrapper provided by numpy_interface.dataset_adapter offers a few convenience functions to traverse composite datasets. Here is a simple example:

>>> for ds in mbw:
>>>    print type(ds)
<class 'vtk.numpy_interface.dataset_adapter.PolyData'>
<class 'vtk.numpy_interface.dataset_adapter.PolyData'>

This wraps up the blog series on numpy_interface. I hope to follow these up with some concrete examples of the module in action and some other useful information on using VTK from Python. Until then, cheers.

Note: This article was originally published on the Kitware blog. Please see the Kitware web site, the VTK web site and the ParaView web site for more information.