Image index ordering in Python API

Engineering
1

I know this topic has been discussed before, and maybe all decisions have been made and too much legacy exists anyway already…

Still, I find it strange to have this numpy array like view of an ITK image with inverted index ordering:

  • coming from C++ and being familiar with the C++ ITK API it seems strange to invert the ordering.
  • I find it error-prone to have an
    • ITK Python API that uses x,y,z-ordering GetPixel(x,y,z), GetDimensions() etc.,
    • a pythonic ITK APi that uses x,y,z-ordering
    • a numpy interface ordering (inspired by scikit-image) with z,y,x ordering
  • it is easy to mix itk.size(img), img.shape, and end up with confusion
  • it gets even more complicated (or not?) when we are dealing with orientation, e.g. img.GetDirection() vs img['direction']

If performance is (was) the main reason to choose this ordering I would like to point out that the analysis linked here, may have not considered wrapping the C++ memory using fortran ordering (basically: reverse strides and shape).

Running the example code from scikit-image, using a 3D image with dimensions 500x500x500, I would apparently get a speed-up of factor 22x for the out of order multiplication. BUT if I do the out of order multiplication on a fortran view of the 3d image, the execution time is nearly the same.

in_order_multiply: 0.12 seconds
out_of_order_multiply: 0.13 seconds
Speedup: 1.1x

Here is the code:

import numpy as np

def test_profile_image_ordering():

    fortran_view = lambda x: x.transpose(np.flip(np.arange(len(x.shape))))

    def in_order_multiply(arr, scalar):
         for plane in list(range(arr.shape[0])):
            arr[plane, :, :] *= scalar

    def out_of_order_multiply(arr, scalar):
        for plane in list(range(arr.shape[2])):
            arr[:, :, plane] *= scalar

    import time
    rng = np.random.default_rng()
    im3d = rng.random((500, 500, 500))
    im3d_f = fortran_view(im3d)
    assert im3d_f.flags.owndata == False

    s0 = time.time(); x = out_of_order_multiply(im3d_f, 5); s1 = time.time()
    print("%.2f seconds" % (s1 - s0))

    t0 = time.time(); x = in_order_multiply(im3d, 5); t1 = time.time()
    print("%.2f seconds" % (t1 - t0))

    print("Speedup: %.1fx" % ((s1 - s0) / (t1 - t0)))

if __name__ == "__main__":
    test_profile_image_ordering()

Maybe the python world prefers this ordering/convention, or maybe it works better with other packages, including some machine learning libraries. I guess there is also the issue of how to create an itk image view of a C-ordered numpy array.
It is a bit unfortunate though for people like me coming from C++ ITK, especially since many simple operations can be implemented very elegantly using numpy array views of the image, and the new pythonic ITK API.

2 Likes
2

I was unhappy about the different conventions for a good while, but David Gobbi’s explanations on the VTK forum helped a lot. The main takeaway for me is this:

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.

2 Likes
3

Thanks @lassoan. Yes, it helps to read that thread. :+1:

4

Just to say that the neuroimaging convention in Python - is to use x, y z indexing. This is the convention in Nibabel, the standard input-output library that nearly every neuroimaging package uses. I wonder whetther the convention in scikit-image starts from the idea of a stack of 2D images, rather than one 3D image.

1 Like
5

It would be nice if the index order was exceedingly apparent at every step. Instead of “image” style vs. “array” style, to me it is “Fortran” index order vs. “C” / “numpy” index order. Perhaps a combination of the two: “ImageF” vs. “ArrayN” in the ITK class, method, and parameter names, and in user variable names in our code samples, would be a constant reminder of which index ordering is appropriate.

“Mommy, why is there an ‘N’ at the end of ‘Array’?” “Well dear, it is to remind us that this object uses the numpy index-ordering convention.” “Oh, I see. Thanks, Mom!”

2 Likes
6

Index ordering conventions are certainly an area of confusion, possible bugs, and performance considerations.

Depending on perspective, many people find one convention or the other strange. Which is strange depends on the perspective. :slight_smile:

There is some documentation here:

https://itkpythonpackage.readthedocs.io/en/master/Quick_start_guide.html

but perhaps we could do better?

We could also encourage the use of xarray.DataArray – conversion is possible with itk.xarray_from_image and itk.image_from_xarray. In this case, axis indexing is explicit with a name, e.g. x. This is more verbose to write but easier to read, reason about, and it avoids errors.

Perhaps we could add itk.size, itk.spacing attribute access via Python @property’s that use NumPy indexing order. Then, to avoid confusion, the rule would be: attribute access, e.g. image.spacing, image.size, image.shape is NumPy/scikit-image index order, function access, e.g. image.GetSpacing() is ITK/Nibabel/VTK index order. What do you think?

numpy.ndarray has a flag that indicates whether a view is Fortran or C order, and numpy operation implementations can support either memory organization. Unfortunately, the assumption is C-ordering in ITK algorithms implementations, and it is not feasible to add support for Fortran ordering.

7

My primary language is C++, not Python, so my opinion might not be widely shared. But this sounds like it would create more confusion than it would solve.

8

While image.shape is in NumPy index order and confusion could be avoided if NumPy changed its convention, this convention cannot be changed. Given this situation, we may have to take the best approach possible, even if it conflicts with what we are used to or how we would like it to be.