Image index ordering in Python API

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.

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.

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

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.

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!”