Why VTK Arrays Are Not NumPy Arrays (Anymore)

In my last blog, I promised to get into the philosophy of the latest rewrite of VTK arrays in Python. Here it is. Let’s start with the line that looks like a step backwards.

>>> import numpy as np
>>> from vtkmodules import vtkImagingCore as imaging_core
>>> w = imaging_core.vtkRTAnalyticSource()()
>>> rtdata = w.point_data['RTData']
>>> isinstance(rtdata, np.ndarray)
False

In 2014, and again in my 2024 take, the same check printed True. Getting a real numpy array was the whole point. So why did we stop?

Three rewrites

The first was numpy_support, which Prabhu Ramachandran added in 2008. It converts: vtk_to_numpy() and numpy_to_vtk(). The second was the numpy_interface module from 2014. Its VTKArray is, to quote its docstring, “a sub-class of numpy ndarray that stores a reference to a vtk array as well as the owning dataset. The numpy array and vtk array should point to the same memory location.”

That last sentence is the catch. It holds as long as the VTK array is one block of memory that numpy can describe. In 2014, that was nearly always true. Today, it often isn’t. Here is an SOA array, the kind we get from simulation codes that store x, y and z separately. I’m using the old dataset_adapter only to show what it does:

>>> from vtkmodules import vtkCommonCore as common_core
>>> from vtkmodules.numpy_interface import dataset_adapter as dsa
>>> x = np.arange(4, dtype=np.float32)
>>> soa = common_core.vtkSOADataArrayTemplate[np.float32]([x, x * 10, x * 100])
>>> old = dsa.vtkDataArrayToVTKArray(soa)
>>> isinstance(old, np.ndarray)
True
>>> old[1, 0] = 42
>>> soa[1, 0]
np.float32(1.0)

The write went nowhere. An ndarray can describe one strided block of memory. This array is three blocks, so old is a copy, and nothing tells you. Ugh.

SOA is the mild case. An implicit array has no memory to point at (an affine array is a slope and an intercept), and an array that lives on a GPU has memory the host can’t touch. For all of these, “a numpy array that points to the same memory” is either a copy or impossible.

Protocols, not inheritance

We subclassed ndarray in 2014 for a reason: it was the only way to get numpy’s behavior. numpy had no way for another object to take part in ufuncs until __array_ufunc__ arrived in NumPy 1.13 (NEP 13), and none for functions such as np.mean() until __array_function__ in NumPy 1.17 (NEP 18). That is also why the 2014 module came with its own algorithms module, with its own max() and friends.

With those protocols in place, we could flip the relationship around. The VTK array is the object. It implements numpy’s protocols, and each array type decides how. The AOS array hands numpy a zero-copy view. The SOA array works one component at a time:

>>> soa[1, 0] = 42
>>> x[1]
np.float32(42.0)
>>> soa * 2
VTKSOAArray(shape=(4, 3), dtype=float32)

The write lands in x, the memory we built the array from, and the result of the arithmetic is still an SOA array. An affine array doesn’t even do the arithmetic. It updates its slope and intercept:

>>> a = common_core.vtkAffineArray[np.float32](1_000_000_000, slope=1)
>>> b = a * 0.5 + 1
>>> b
VTKAffineArray(slope=0.5, intercept=1.0, num_values=1000000000, dtype=<class 'numpy.float32'>)
>>> b.GetActualMemorySize()
1

A billion values in 1 KiB. Try that with an ndarray. Composite datasets get the same treatment. One array spans all the blocks without concatenating them:

>>> from vtkmodules import vtkCommonDataModel as data_model
>>> pds = data_model.vtkPartitionedDataSet()
>>> pds.append(w)
>>> pds.append(imaging_core.vtkRTAnalyticSource(whole_extent=(0, 4, 0, 4, 0, 4))())
>>> pds.point_data['RTData']
VTKPartitionedArray(shape=(9386,), dtype=float32, blocks=2)

As an ndarray, that would have to be a copy too.

Arrays that know where they came from

There is a second reason, and I wrote about it back in 2014: “VTK brings a full-featured data model that goes beyond n-dimensional arrays to the picture.” An array from a dataset knows which dataset it came from, and so does the result of a mesh-aware function:

>>> rtdata.dataset is w
True
>>> from vtkmodules.util import functions
>>> grad = functions.gradient(rtdata)
>>> grad.shape
(9261, 3)
>>> grad.dataset is w
True

gradient() needs the mesh, and it finds it through the array. The 2014 VTKArray carried the same information, but only if the data went through WrapDataObject() first. Now every array you get from point_data has it. There is a lot more to say about functions. Another blog.

The price

Nothing is free. Here is what this costs.

• Code that checks isinstance(x, np.ndarray) says no. For those, to_numpy() gives you a real ndarray: a view for AOS arrays, a copy for SOA and implicit arrays, which have no single block of memory to share.
• We now own numpy’s semantics. When you implement __array_ufunc__, you implement numpy’s rules, fine print included. While writing my last blog, I found that rtdata + 2 came out float64 where numpy keeps float32. The fix then failed on Windows, because VTK’s Windows CI runs numpy 1.26, whose rules differ from numpy 2’s. Both are fixed now (VTK 9.7.2).
• Lazy isn’t bit-identical. An affine array computes its new slope and intercept in double precision, so a float32 result can differ from numpy’s in the last digit or so.

I think this is the right trade. The alternative is an ndarray that is silently a copy in exactly the cases VTK cares about.

Where this goes

Protocols over inheritance is also where the rest of the Python world went. The Python array API standard is built around __array_namespace__(), which VTK arrays don’t implement yet. GPU libraries such as PyTorch and CuPy share memory through DLPack, and VTK has started to speak it too. Both fit the “VTK array is the object” design much better than a subclass ever could.

In my next blog, I will look at something I glossed over in “The price”: when VTK and numpy share memory, who owns it, and what happens when one side lets go?

Many thanks to Prabhu Ramachandran for the initial implementation of numpy_support in 2008.

Are you working with VTK in your scientific computing applications? Kitware’s experts can help you address technical challenges and develop solutions tailored to your needs. Contact us to discuss your project.