Why to invert the order of the indices between numpy and SimpleITK?

[Diego_C (Diego Cantor)]

HI,

Can anyone explain the logic between this:

This is taken from the SimpleITK Notebooks

It states that the order of the indices is the opposite in numpy. But in fact this is an arbitrary choice. For instance in matplotlib the last index of the numpy array represents the depth:

This type of change back and forth between channels first and channels last is prone to errors if you forget to switch the order of your indices and also reduces interoperability with other libraries (like matplotlib).

Could I suggest to have a way to configure your preferred order for numpy arrays in a similar fashion to other libraries such as keras:

https://keras.io/backend/

[blowekamp (Bradley Lowekamp)]

So imshow does correctly work with the layout that SimpleITK exports. Here is an RGB example:

In [14]: img = sitk.ReadImage("simpleitk_logo.png")

In [15]: nda = sitk.GetArrayViewFromImage(img)

In [16]: nda.shape
Out[16]: (53, 200, 4)

In [17]: img.GetSize()
Out[17]: (200, 53)

In [18]: plt.imshow(img)

The RGB channels ARE last.

The notebook example may be a little confusing because it is utilizing the size 3 for the z-physical space dimension, it is not an RGB image e.g. a 2-d image with 3 components. Additionally, imshow does not support 3-D spacial images, only 2-D image with 1,3, or 4 channels. The term depth in ambiguous here either referring to the number of channels or the 3rd spacial dimension.

The choices were not arbitrary. SimpleITK was designed to be consistent with ITK. ITK’s indexes are represented as x,y,z, so SimpleITK must reference the images that way too. The GetArrayViewFromImage method does not copy the ITK image buffer, so the layout of the data does not change when exporting. Numpy array created is index via C conventions (buffer[z][y][x][c]) vai row major with the fastest axis being the last “index”. So this approach is both efficient and compatible with numpy.

But you are correct it is error prone and confusing.

[Diego_C (Diego Cantor)]

I don’t disagree with you for a second that SimpleITK needs to be compatible with ITK. Why in heavens would not be? That wasn’t the spirit of my question.

I am thinking about interoperability between SimpleITK and Numpy.
What I am thinking is why do we need to flip the order of the indices when getting the numpy arrays??

Notice that in your example you pass the sitk image to imshow, it makes sense that it works since the internal layout is [x,y,z] as show consistently throughout the sitk image methods. So no problem there.

But there is the expectation operate in the same manner with the numpy array, but you can’t because the indices are reversed.

So if I were going to extract the numpy array and do any kind of operation and then visualize it, you have to be constantly aware of this change or your numpy processing code won’t work or neither the visualization.

Maybe this [z, y, x] ordering is compatible with numpy as you suggest but so are other orderings like [z,x,y] or [y,z,x] but I would debate that computational efficiency is relative to the way the programmer handles his/her arrays.

Every time I pass from sitk to numpy and back I need to be aware of this change. That reduces compatibility, and creates the perfect situation to introduce bugs in your code.

My suggestion is to have an API method that enables you to return numpy arrays with the order that you expect, in a similar fashion to the example I provided earlier for keras.

Something like:

...
sitk.SetArrayLayout(sitk.DEPTH_FIRST)
narr = sitk.GetArrayFromImage(img) # returns array as [z,y,x]

sitk.SetArrayLayout(sitk.DEPTH_LAST)
narr = sitk.GetArrayFromImage(img) # returns array as [x,y,z]

[phcerdan (Pablo Hernandez-Cerdan)]

No idea about the python bindings so I cannot comment much on this, but if GetArrayFromImage is not performing a copy as @blowekamp pointed, I hardly see how this can be implemented. Maybe a method doing a copy with whatever layout CopyArrayFromImage(img, layout), as you suggested can be introduced, but it won’t be able to inter-operate between numpy and itk without an extra copy back to ITK, CopyImageFromArray(array, layout) which would be inefficient, but might still be useful.

@Diego_C

I don’t disagree with you for a second that SimpleITK needs to be compatible with ITK. Why in heavens would not be? That wasn’t the spirit of my question.

Speaking as a fellow user/consumer of ITK, I recommend you to avoid this kind of comments. We all want the toolkit to be the sharpest tool possible, and we do it in a collective manner, in our time, with step-by-step improvements. Stick to pro-active, helpful comment and the continuous process of improvement will be easier and more enjoyable for everybody.

[phcerdan (Pablo Hernandez-Cerdan)]

@phcerdan

CopyImageFromArray(array, layout)

This can be improved with CopyImageFromArray(array, layout, itk_image_metadata)
To create the image with the metadata (origin, spacing, direction) from an existing itk image.

And/or OverrideImageFromArray(array, layout, image) to change the data of the image in-place.

[lassoan (Andras Lasso)]

@blowekamp

Numpy array created is index via C conventions (buffer[z][y][c]) vai row major with the fastest axis being the last “index”. So this approach is both efficient and compatible with numpy

You can specify C/Fortran convention for numpy.asarray and reshape. Could those be used to get consistent indexes?

[Diego_C (Diego Cantor)]

@phcerdan

Maybe a method doing a copy with whatever layout CopyArrayFromImage(img, layout), as you suggested can be introduced, but it won’t be able to inter-operate between numpy and itk without an extra copy back to ITK, CopyImageFromArray(array, layout) which would be inefficient, but might still be useful.

Hi Pablo!!

I thought of this as well, but the problem to have another method to do something that is already been done by GetArrayFromImage is that now the programmer has two (well, four!) methods to learn and keep track of. So I think a global method for the programmer to decide what order he/she expects from the Numpy arrays returned by SimpleITK (similar to the channels-first, channels-last options in keras) would be useful:

Set it once, at the beginning of your code, so you don’t have to deal with this decision every time you get a numpy array from SimpleITK.

Question: if GetArrayFromImage is not making a copy… Does that mean that any changes in the data contained by the numpy array, will also be shared if I query the image directly with GetPixel or similar methods?

Thanks,

Diego

[blowekamp (Bradley Lowekamp)]

@Diego_C

Question: if GetArrayFromImage is not making a copy… Does that mean that any changes in the data contained by the numpy array, will also be shared if I query the image directly with GetPixel or similar methods?

GetArrayFromImage makes a deep copy, but GetArrayViewFromImage is a read only buffer in numpy. With the view you must maintain a ssitk::Image copy and if you modify the image it could effect the bumpy buffer. There is also SimpleITK’s lazy copying and reference accounting which can affect which image the buffer points to.

[lassoan (Andras Lasso)]

For me it seems that we should be able to get IJK-indexed view (see my comment above, using Fortran ordering) and then everything would work without manually reversing axes.
@blowekamp Have you tried to use Fortran ordering?

[fbudin (Francois Budin)]

Based on the little I remember when I implemented this and what I can see in the implementation of _GetImageViewFromArray() which allows to import images in ITK from NumPy arrays that have Fortran or C ordering, it should be possible to add an option to _GetArrayViewFromImage() to select the ordering in the imported NumPy array. Does anybody want to give it a shot? Potentially it is just a matter of changing the order of the size vector components.

[blowekamp (Bradley Lowekamp)]

The current approach is correct. It follows ITK conventions ( and ITK Python ), and numpy conventions. ITK index via [i,j,k] and numpy index via [k,j,i]. Look at other Python imaging libraries such as PIL or scikit-image.

Yes you can utilize numpy.reshape to convert to fortran order while reversing the shape array, but that may result in a copy. And when you apply algorithms to the fortran array it may have to copy the data so that the buffer is in the proper order. This is not recommended.

Regarding keras, you have conflated index order and image data layout options. keras is following numpy conventions for [k,j,i] indexing. For RGB ( multi-components ), it has the option for the RGB to be the fast axis or the slow axis. For the former you can consider the layout as each pixel is a RGB value, while with the latter the data is laid out as the R-image, then G-image, then B-image.

[lassoan (Andras Lasso)]

Thanks for the explanation, it makes sense. We want to keep the common convention of indexing matrices as (row, column), but unfortunately for images it means (y, x) coordinates - kind of inverted. There is nothing else we can do about this, people just have to be aware and pay attention. From the variable type (image or matrix) it is always clear what indexing convention is used.

[Diego_C (Diego Cantor)]

@blowekamp

GetArrayFromImage makes a deep copy, but GetArrayViewFromImage is a read only buffer in numpy.

HI Bradley, thanks for the explanation. I was referring indeed to GetArrayFromImage not to GetArrayViewFromImage, so if it does a deep copy as you say, it would be possible then to return that copy in the expected [x,y,z] (as in SimpleITK and ITK) order, is that correct?

Thanks,

Diego

[Diego_C (Diego Cantor)]

I checked on PIL and scikit-image as per @blowekamp suggestion. Bradley is correct in that scikit-image handles operations on 3D images using C-contiguous arrays, where as channels in RGB images go last.

In case you guys are interested, this link explains very clearly the difference between RGB and 3D images in terms of the array layout for scikit-image 4. A crash course on NumPy for images — skimage 0.23.0rc0.dev0 documentation

An excerpt explaining everything is here:

Perhaps adding some like this to SimpleITK’s documentation could be really helpful? What do you guys think?

On the other hand I didn’t find info on PIL other than when you retrieve data from the image, you get a cordinate-less array in the sense that it is a flattened array instead of multi-indexed.

Regards,

Diego

[zivy (Ziv Yaniv)]

@Diego_C, we intend to work on general improvements of the SimpleITK documentation in the next couple of months. Thanks for the pointer to the scikit-image documentation, this is useful for us as a point of reference.

[Jingnan_Jia]

Hi, Diego_C, I noted that in scikit, the order of coordinates is (row, col) or (x, y) no matter for 2D or 3D images. However, in SimpleITK, the order changed from (x, y) to (y, x) after GetArrayFromImage. Could you explain if it is correct? Why x and y need to be exchanged?