Numpy: Striding a multiple channel image

Question

I'm having some trouble getting my head around striding in numpy. I'm writing some code that does interpolation of multi-channel images. I define my images as 3 dimensional arrays of type np.ndarray with shape [HEIGHT x WIDTH x CHANNELS]. The C++ that I am writing must work in both Matlab AND Python. For single channel images, my code works fine, and for multi-channel images in Matlab, my code works fine.

In order to interpolate an image, I am writing a method whereby, given an [M x N x P] array, you can provide a set of X and Y sub-pixel coordinates to be interpolated within the image. This is identical to the functionality of scipy's ndimage.map_coordinates. Unfortunately, I require an interpolation method that yields identical results in both Matlab and Python and am thus rolling my own interpolation code.

My issue is that Matlab arranges it's 3-dimensional memory by stacking the concatenating the channels one after another. This means that, for a [10, 10, 2] image, the first 100 elements will be the first channel, and elements [100, 200] elements will be the second channel. Therefore, to index in to the Matlab contiguous memory, I index as follows:

// i is the element of the indices array
// j is the current channel
// F is the image we are indexing
// F_MAX is M * N (the number of pixels per channel)
// N_ELEMS is the total number of elements in the indices array
// f_index is the index in the contiguous array equivalent to the x and y coordinate in the 2D image
for (size_t j = 0; j < N_CHANNELS; j++)
{
  out[i + j * N_ELEMS] = F[f_index + j * F_MAX];
}

My issue is that numpy orders it's 3-dimensional arrays along the 3rd axis. That is to say, given a [10, 10, 2] array, the first 2 elements are indices [0, 0, 0] and [0, 0, 1]. In Matlab they are indices [0, 0, 0] and [0, 1, 0].

I think I can rectify my issue by using a stride in numpy. However, I am totally failing to come up with an appropriate stride pattern. Therefore, for my example of a [10, 10, 2] array, how can I change the strides, from (assuming doubles):

>>> np.ones([10,10,2], dtype=np.float64).strides
(160, 16, 8)

to something that I can index into as I do for Matlab arrays?

I should mention that I am aware of the column major/row major difference between Matlab and numpy respectively. As stated, my method works for single channel images but indexes wrong with more than 1 channel.

Answer 1

Perhaps you can make use of the function np.swapaxes as in the following ipython example:

In [1]: a = np.arange(2*2*2).reshape((2,2,2))

In [2]: a
Out[2]: 
array([[[0, 1],
        [2, 3]],

       [[4, 5],
        [6, 7]]])

In [3]: a.flatten()
Out[3]: array([0, 1, 2, 3, 4, 5, 6, 7])

In [4]: np.swapaxes(a,2,1).flatten()
Out[4]: array([0, 2, 1, 3, 4, 6, 5, 7])

EDIT:

I think the internal memory layout only changes after you take a copy of the array with swapped axes, see:

In [6]: b = a.swapaxes(1,2)

In [7]: b.strides
Out[7]: (16, 4, 8)

In [8]: b = a.swapaxes(1,2).copy()

In [9]: b.strides
Out[9]: (16, 8, 4)