Who Owns This Memory?

In my last blog, I said that to_numpy() gives you a view for some VTK arrays and a copy for others. That deserves
a closer look. When two libraries share memory, the interesting question is always the same: who owns it? Who
frees it, and when? Let’s start on the numpy side.

From numpy to VTK

>>> import gc
>>> import numpy as np
>>> from vtkmodules import vtkCommonCore as common_core
>>> a = np.arange(6, dtype=np.float64)
>>> v = common_core.vtkDoubleArray(a)
>>> a[0] = 100
>>> v[0]
np.float64(100.0)
>>> np.shares_memory(v, a)
True

The VTK array doesn’t copy a. It points at numpy’s memory, so a change on one side shows up on the other. Now
the fun part:

>>> del a
>>> gc.collect() >= 0
True
>>> v[0]
np.float64(100.0)

Still there. numpy owns the memory, but the VTK array’s buffer holds a reference to the numpy array, so numpy
can’t free it as long as v is alive. You don’t have to keep a around yourself.

Datasets work the same way, with one catch:

>>> from vtkmodules import vtkImagingCore as imaging_core
>>> w = imaging_core.vtkRTAnalyticSource()()
>>> b = np.zeros(w.number_of_points)
>>> w.point_data['zeros'] = b
>>> np.shares_memory(w.point_data['zeros'], b)
True

>>> c = np.zeros((w.number_of_points, 3))[:, 0]
>>> w.point_data['strided'] = c
>>> np.shares_memory(w.point_data['strided'], c)
False
>>> c[0] = 1
>>> w.point_data['strided'][0]
np.float64(0.0)

c is a column of a 2D array, so its values are not next to each other in memory. A VTK array needs them to be
(well, most of them do, more on that below), so the assignment makes a contiguous copy. Nothing tells you, and
changes to c no longer reach the dataset. If sharing matters, check with np.shares_memory() or hand over a
contiguous array.

Assigning a VTK array is a slightly different story. As I showed before, the dataset gets a new VTK array object that
points to the same memory. Same memory, different object.

From VTK to numpy

Going the other way, it depends on the array. For a regular (AOS) array, to_numpy() is a zero-copy view, and you
can write through it:

>>> rtdata = w.point_data['RTData']
>>> view = rtdata.to_numpy()
>>> np.shares_memory(view, rtdata)
True
>>> view.flags['WRITEABLE']
True

An SOA array keeps each component in its own block of memory, and a single ndarray can’t describe three blocks.
So to_numpy() stitches them together into a copy. If you want zero copy, ask for one component at a time:

>>> x = np.arange(4, dtype=np.float32)
>>> soa = common_core.vtkSOADataArrayTemplate[np.float32]([x, x * 10, x * 100])
>>> np.shares_memory(soa.to_numpy(), x)
False
>>> comp0 = np.asarray(soa.GetComponentBuffer(0))
>>> np.shares_memory(comp0, x)
True

Implicit arrays have no memory at all, so to_numpy() computes the values into a fresh array. Writing to that array
does nothing to the VTK array.

>>> aff = common_core.vtkAffineArray[np.float64](4, slope=2)
>>> aff.to_numpy()
array([0., 2., 4., 6.])

Back in 2015, VTK at least printed a warning when GetVoidPointer() had to generate a copy of a mapped array.
to_numpy() copies an SOA array without a word. That is fine, as long as you are aware.

Who keeps it alive

In VTK, the object that owns an array’s memory is a vtkBuffer, and as of VTK 9.7 it supports Python’s buffer
protocol. That makes it the lowest level you can go:

>>> buf = rtdata.GetBuffer()
>>> raw = np.asarray(buf)
>>> raw.shape, raw.dtype
((9261,), dtype('float32'))
>>> np.shares_memory(raw, rtdata)
True
>>> mv = memoryview(buf)
>>> mv.format, mv.itemsize, mv.nbytes
('f', 4, 37044)

The buffer protocol is also what makes views safe to keep around. Let’s throw away the dataset and the algorithm
that made it, and hold on to a view:

>>> src = imaging_core.vtkRTAnalyticSource()
>>> image = src()
>>> keep = image.point_data['RTData'].to_numpy()
>>> first = float(keep[0])
>>> del image, src
>>> gc.collect() >= 0
True
>>> float(keep[0]) == first
True

The view is fine. It holds on to the vtkBuffer (through keep.base, a memoryview), and the vtkBuffer is reference
counted, so the memory lives as long as the view does. The dataset can go away. The memory can’t.

When memory moves

There is one case where this breaks down. A VTK array can grow, and when it outgrows its memory, it reallocates.
Sometimes the block grows in place and nothing bad happens. Sometimes it moves: VTK gets a new block, copies
the values over and frees the old one. Any numpy view is then left pointing at freed memory.

>>> d = common_core.vtkDoubleArray()
>>> d.SetNumberOfValues(3)
True
>>> view = d.to_numpy()
>>> view[:] = [0, 1, 2]
>>> d.InsertNextValue(3.0)
3

and you get this warning:

RuntimeWarning: The underlying VTK array has reallocated its buffer. Any numpy
view previously obtained from this array is now stale and points to invalid
memory. Please retrieve a fresh view from the array.
>>> np.shares_memory(view, d)
False
>>> view = d.to_numpy()
>>> view
array([0., 1., 2., 3.])

This time the block moved. VTK arrays fire a new event (new in VTK 9.7), BufferChangedEvent, whenever they
reallocate, and the Python layer listens for it. That’s where the warning comes from. From Python, you can’t tell
whether the block moved or grew in place, so the warning fires every time, and it’s a warning, not a fix. The stale
view is still yours, and reading it after VTK has freed that memory is undefined behavior. So the rule is simple: if
you resize an array (InsertNextValue(), SetNumberOfTuples(), Resize() etc.), get a fresh view afterwards.

In my next blog, I will look at arrays that sidestep the question altogether, because they have no memory at all:
implicit arrays.