Very confused about imdata matrix index order

Support
1

Dear,
I am working in python. I am very confused about the index ordering in the matrix I do pass to imdata.
I do have a numpy array named matrix, with indexes [x,y,z] that I pass to a vtkImageData with:

imdata = vtk.vtkImageData()
depthArray = numpy_support.numpy_to_vtk(
         matrix.ravel(), deep=True,
         array_type=numpy_support.get_vtk_array_type(matriz.dtype)
     )
imdata.SetDimensions( matrix.shape[0],matrix.shape[1],matrix.shape[2])
imdata.SetSpacing(1,1,1)
imdata.GetPointData().SetScalars(depthArray)

but when I plot it, it looks rotated in space. I guess the ordering [x,y,z] in the matrix is not correct. Is it ok?

Is x,y,z the correct index order for other functions?

Regards,
Daniel

1 Like
2

For an array (or matrix), the indexing is [row, col] for two dimensions, or [slice, row, col] for three dimensions.

Consider the following code, where a[1,0] means row=1, col=0:

>>> import numpy as np
>>> np.array([[1,2],[3,4]])
array([[1, 2],
       [3, 4]])
>>> a = np.array([[1,2],[3,4]])
>>> a[1,0]
3

Since arrays are indexed as [slice,row,col] and slice=z, row=y, col=x it is natural that, if the array stores an image, we index the array as follows:

>>> a[Z, Y, X] # 3D
>>> b[Y, X] # 2D

Here is the best way to think about it: the order of the indexes is according to the type of object the data is stored in. A numpy array has an ‘array-style’ interface, and is indexed by [row, col] or [slice, row, col]. A vtkImageData has an ‘image-style’ interface, and is indexed by (x, y) or (x, y, z). So that is how VTK works.

You have to be extra-careful when using python image libraries like nibabel, because they transpose the array in order to use image-style indexing with the arrays. In other words, they use ‘row=x’ and "col=y’ which seems unnatural to me.

3 Likes
3

Thanks!

And the order in commands is like:
VtkImageData.SetSpacing(x,y,z)
?

I guess I do have a mess with numpy.meshgrid which I use to generate values for the matrix calculations.

Regards,
Daniel

4

Yes.

The vtkImageData index order is (x,y,z).
The numpy array index order is [row, col] or [slice, row, col].

So when data is moved from a numpy array into a vtkImageData, the index order might change.
The indexing order is determined by the container’s interface.

1 Like
5

Thanks again!

6

In numpy’s terms, VTK expects Fortran-order, while numpy defaults to C-order for its memory layout.

I typically use np.meshgrid with indexing='ij', and pass order='F' when reshaping the data to VTK format (matrix.ravel(order='F') in your example).
A full example can be found here:

1 Like
7

VTK (and ITK) always use C-contiguous memory layout. I think you are confusing index ordering and memory ordering.

Calling matrix.ravel(order=‘F’) on a C-contiguous numpy array will generate a transposed copy of the data. It’s an inefficient way of dealing with the difference in index ordering, so I would not recommend it.

As I tried to explain in my earlier posts, the issue is differences in interface (not memory order). Since vtkImageData was created for people who work with 3D graphics, the most natural index ordering is (image column, image row, image slice). Numpy was originally created for linear algebra, so its interface uses [matrix row, matrix column] with the possibility of additional indices either before or after. So I’ll repeat myself like a broken record: the difference in index order between VTK and numpy is due to a difference in interface rather than a difference in the memory order of the data.

I think a lot of the confusion arises because numpy is very flexible and very abstract. When you store an image in a numpy array, you get to choose whether the numpy indices correspond to [x,y,z] or [z,y,x] or even [y,z,x]. And you can also, separately, choose the memory order by choosing the strides. VTK offers no such flexibility, since vtkImageData is always (x,y,z) and always C-contiguous.

If you really want to use [x,y,z] ordering with numpy and smoothly transition your numpy arrays into vtkImageData without doing extra copies of the data, your best bet is to keep all your numpy arrays F-contiguous. My own preference, however, is to is to use C-contiguous numpy arrays and index them with [z,y,x].

Edit: On ITK discourse I saw a useful link about coordinate conventions in scikit-image. The recommendations on that page are also useful for people who want to use numpy and VTK.

2 Likes
9

Thanks! My problem was the order in numpy.meshgrid, I was generating the image in a different orientation that the needed

10

Thanks @dgobbi your explanations are really helpful (even rephrasing them in slightly different ways adds value).

However, I’m not sure what you mean by your comment about VTK having no flexibility in memory layout. Wouldn’t the "new" array dispatch mechanism allow alternative memory layouts? Or you mean that that’s a theoretical possibility because in practice only one layout is supported?

11

Hi Andras, the VTK image filters cannot use the VTK SOA arrays, since they require a raw memory pointer to a contiguous block of memory. It’s similar to the situation with image orientation: even though it’s a core feature of VTK, most VTK image filters do not implement it, and as far as I know it isn’t on any roadmap. It’s a topic worthy of discussion, but in a separate thread.

AOS versus SOA is not quite the same as C-contiguous versus F-contiguous, since SOA arrays do not have to be contiguous. Although it is theoretically (and practically) possible to handle F-contiguous data via SOA arrays.

Edit: Also, vtkImageData interprets the underlying array dimensions as (tuple index, component index) rather than as spatial dimensions. So SOA vs AOS in images would correspond to ‘planar’ vs. ‘packed’.

12

Thanks for the explanation. I could not understand all the details (e.g., by “SOA vs AOS in images would correspond to ‘planar’ vs. ‘packed’” do you mean IJKC vs CIJK?), but overall it seems that although the array dispatch mechanism provides some flexibility, in practice only a single layout works, because alternatives are not implemented.

This is in 3D Slicer’s roadmap. We have a couple of classes (such as vtkBSplineTransform, vtkGridTransform, vtkPolyDataToImageStencil) extended to support image orientation, but did not have time to integrate it into VTK yet. A complete transition to use image orientation stored in vtkImageData in Slicer and make most commonly used VTK classes to use image orientation should happen in the next few years.

1 Like