Arrays That Aren’t There
In my last blog, I said that to_numpy() on an implicit array computes the values into a fresh array, because an implicit array has no memory. It’s time to look at them properly. You have probably used one already without knowing it:
>>> import numpy as np
>>> from vtkmodules import vtkImagingCore as imaging_core
>>> from vtkmodules import vtkCommonCore as common_core
>>> w = imaging_core.vtkRTAnalyticSource()()
>>> pts = w.points.data
>>> pts
VTKStructuredPointArray(shape=(9261, 3), dims=(21, 21, 21), dtype=float64)
>>> pts.GetActualMemorySize()
1
>>> pts[0], pts[-1]
(array([-10., -10., -10.]), array([10., 10., 10.]))9261 points with 3 doubles each would take about 217 KiB. This array takes 1 (GetActualMemorySize() reports KiB). An image knows its origin, spacing and dimensions, and the point array computes a coordinate when you ask for it.
Back in 2015, I wrote a constant array and an implicit points array for image data as examples of what VTK’s then-new mapped arrays could do. Both ideas are in VTK proper now, as vtkConstantArray and vtkStructuredPointArray, along with a few more, and they all behave like numpy arrays in Python.
The family
Some implicit arrays compute their values from a formula:
>>> a = common_core.vtkAffineArray[np.float64](1_000_000_000, slope=0.001)
>>> a
VTKAffineArray(slope=0.001, intercept=0.0, num_values=1000000000, dtype=<class 'numpy.float64'>)
>>> a.GetActualMemorySize()
1
>>> c = common_core.vtkConstantArray[np.float64]((1000000, 3), 1.5)
>>> c
VTKConstantArray(value=1.5, shape=(1000000, 3), dtype=float64)A billion values (think numpy.arange()) in 1 KiB, and a million 3-component tuples of the same value, also in 1 KiB. Others are views of arrays that do have memory. An indexed array picks values out of another array, in any order, without copying them:
>>> base = common_core.vtkFloatArray(np.array([10, 20, 30, 40, 50], dtype=np.float32))
>>> ids = common_core.vtkIdTypeArray(np.array([4, 0, 2]))
>>> sub = common_core.vtkIndexedArray[np.float32](ids, base)
>>> sub.to_numpy()
array([50., 10., 30.], dtype=float32)
>>> base.SetValue(0, 99)
>>> sub.to_numpy()
array([50., 99., 30.], dtype=float32)A composite array glues arrays together, again without copying:
>>> first = common_core.vtkFloatArray(np.arange(3, dtype=np.float32))
>>> second = common_core.vtkFloatArray(np.arange(3, 6, dtype=np.float32))
>>> both = common_core.vtkCompositeArray[np.float32]([first, second])
>>> both.to_numpy()
array([0., 1., 2., 3., 4., 5.], dtype=float32)And a strided array takes every n-th value of a buffer:
>>> interleaved = common_core.vtkFloatArray(np.arange(12, dtype=np.float32))
>>> every_third = common_core.vtkStridedArray[np.float32](interleaved.GetBuffer(), 3, 0)
>>> every_third.to_numpy()
array([0., 3., 6., 9.], dtype=float32)
The indexed one is the most useful of the views in my experience. It’s how you hand a filter a subset of an array, or the same array in a different order, without paying for a copy.
Staying lazy
Here is what makes them pleasant to use from Python: when there’s a shortcut, the result stays implicit.
>>> b = a * 2 + 5
>>> b
VTKAffineArray(slope=0.002, intercept=5.0, num_values=1000000000, dtype=<class 'numpy.float64'>)
>>> a[::1000]
VTKAffineArray(slope=1.0, intercept=0.0, num_values=1000000, dtype=<class 'numpy.float64'>)
>>> a.sum(), a.max()
(np.float64(499999999500000.0), np.float64(999999.9990000001))
>>> c * 2
VTKConstantArray(value=3.0, shape=(1000000, 3), dtype=float64)
>>> c.mean(axis=0)
array([1.5, 1.5, 1.5])Arithmetic with a scalar gives a new slope and intercept, slicing gives a shorter affine array, and the sum of a billion values comes from a formula instead of a loop. When there is no shortcut, you get a regular array:
>>> np.sin(common_core.vtkAffineArray[np.float64](4, slope=1))
VTKAOSArray(shape=(4,), dtype=float64)(Please don’t try np.sin(a) on the billion-value one. That is 8 GB.)
Where they come from
VTK makes some of them for you, like the image points above. You can make them yourself, as we did. And you can ask VTK to find them for you: vtkToImplicitArrayFilter looks at the arrays you select and replaces the ones it can with implicit arrays.
>>> from vtkmodules import vtkFiltersReduction as reduction
>>> w.point_data['ramp'] = np.arange(w.number_of_points) * 0.5
>>> w.point_data['ramp'].GetActualMemorySize()
73
>>> squeeze = reduction.vtkToImplicitArrayFilter(strategy=reduction.vtkToAffineArrayStrategy())
>>> squeeze.point_data_array_selection['ramp'] = True
>>> squeeze.point_data_array_selection['RTData'] = True
>>> out = squeeze(w)
>>> out.point_data['ramp']
VTKAffineArray(slope=0.5, intercept=0.0, num_values=9261, dtype=<class 'numpy.float64'>)
>>> out.point_data['ramp'].GetActualMemorySize()
1
>>> type(out.point_data['RTData']).__name__
'VTKAOSArray_vtkTypeFloat32Array'The ramp became an affine array, and RTData, which is not a straight line, was left alone. There are strategies for constant arrays and more.
The fine print
Every access is computed. Julien Fausty, Charles Gueunet, Timothée Chabat and Mathieu Westphal measured the cost in C++: close to nothing through GetValue(), 15-20% through VTK’s array ranges, 3-4x for a composite of two arrays, and up to 5x for an indexed array read in random order. Back in 2015, I promised to dig into the cost of these abstractions and never got around to it. They did.
They are read-only.
>>> a[0] = 1
Traceback (most recent call last):
...
TypeError: VTKAffineArray is read-onlyLazy isn’t bit-identical. An affine array computes its new slope and intercept in double precision. For float32 arrays, that can differ from what numpy gets element by element in the last digit.
nbytes is what the values would take. It follows numpy’s meaning. GetActualMemorySize() is what they do take:
>>> c.nbytes, c.GetActualMemorySize()
(24000000, 1)In my next blog, I will get to something numpy cannot do at all: math that needs the mesh, like gradients, vorticity and cell volumes.
Many thanks to Julien Fausty, Charles Gueunet, Timothée Chabat and Mathieu Westphal for the C++ implicit array framework that all of this is built on.