 # Very confused about imdata matrix index order

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,matrix.shape,matrix.shape)
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

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.

2 Likes

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

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

Thanks again!

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

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

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