Multidimensional array shape in Matlab

[naegelejd]

TLDR:
Given the following example yardl types, what is the desired shape of each corresponding multidmensional array in Matlab?

ndarray: double[channel, z, y, x]

matrix: double[row, col]

arrayOfMatrices: double[count, row, col]

---

Yardl currently serializes and deserializes multidimensional arrays in C-order (the last index varies the fastest as you traverse the array).
Matlab, of course, uses Fortran order (the first index varies the fastest as you traverse the array).

This means that the most efficient way for yardl to deserialize a multidimensional array in Matlab is to first reverse the dimensions as specified in the model. For example, given the model:

ImageData: double[channel, z, y, x]

or, using yardl's expanded syntax:

ImageData: !array
  items: double
  dimensions:
    channel:
    z:
    y:
    x:

In Python and C++, elements of an ImageData data are accessed using data[channel, z, y, x].
In Matlab, however, elements of data are accessed using data[x, y, z, channel].
This reversed order of dimensions is "correct" in terms of indexing multidimensional arrays in Matlab.
Thus, it is reasonable to establish a convention in yardl that Matlab N-dimensional arrays always have reversed dimensions with respect to the model.

However, this may lead to user errors when operating on matrices.
A matrix always has two dimensions: rows and columns.
In C++, Python, and Matlab, the entry in the ith row and jth column is at matrix[i, j]

Example:

matrix[r:3, c:5] has 3 rows and 5 columns...

0   1   2   3   4 
5   6   7   8   9
10  11  12  13  14

If we define this matrix in Python and serialize it:

matrix = np.array(range(15), dtype=float).reshape((3, 5))
with open(filename, "wb") as f:
    f.write(matrix.data)

Then, to get the same matrix in Matlab, we need to read the array in reverse dimension order and transpose it:

f = fopen(filename, "r");
d = fread(f, [5 3], 'double');
m = transpose(d)
fclose(f);

    m =

            0     1     2     3     4
            5     6     7     8     9
            10    11    12    13    14

Transposing the matrix in Matlab breaks the "convention" that multidimensional arrays have reversed dimensions.

Likewise, if a user serializes an array of matrices:

ArrayOfMatrices: double[count, row, col]

• Following the convention, in Matlab an ArrayOfMatrices data would have the shape [col, row, count]. The inner matrices may need to be transposed prior to any matrix operations.
• If yardl transposes the full ndarray (permute with dimensions reversed again), it would have the shape [count, row, col], but this is counterintuitive in Matlab, since the count dimension (which represents pages of matrices) should be the last dimension.
• Yardl could permute the ndarray so that its shape becomes [row, col, count], but this would result in the original example ImageData ndarray having shape [y, x, z, channel] which, again, seems confusing.

So, given the following example yardl types, what is the desired shape of each corresponding multidimensional array in Matlab?

ndarray: double[channel, z, y, x]

matrix: double[row, col]

arrayOfMatrices: double[count, row, col]

[ojosephs]

My initial feeling is that, for Matlab, "seems confusing" seems the most appropriate option(!)

So, for a Matlab variable named "data":

ndarray: double[channel, z, y, x]
data( y, x, z, channel)

matrix: double[row, col]
data( row, col)

arrayOfMatrices: double[count, row, col]
data( row, col, count)

My reasoning is that it is very common, in Matlab, for the first index to refer to a "vertical" coordinate (e.g. "row" or "y") and the second index to be the "horizontal" coordinate (e.g. "column" or "x"). This holds for matrices, matrix pages, 2D images and 3D volumes formed from arrays of 2D images. The following help pages have some relevant examples.

https://uk.mathworks.com/help/matlab/ref/pagetranspose.html
https://uk.mathworks.com/help/images/ref/montage.html
https://uk.mathworks.com/help/medical-imaging/ref/medicalimage.montage.html
https://uk.mathworks.com/help/matlab/ref/image.html

[roopchansinghv]

I would lean towards not breaking what users currently expect to get from the ISMRMRD's MATLAB API.

I feel like debugging an array dimension order bug is likely to cause more headaches ...

If there are compelling performance reasons to change this (i.e. the expected output of the MATLAB API), then that change would have to be pretty prominently high-lighted ...

My (inflation-adjusted) $0.03 ...

[dchansen]

For the original ISMRMRD, I remember we did try to keep it consistent between Matlab and Python - at that time, by making python use the Fortran order. We later came to regret that decision, as it just made "Pythonic" code really hard to write.
I think the same holds for Matlab. I would stick with the transposed order.

[rajramasawmy]

Would this convention break existing gadgetron-matlab reading? I think currently the data comes in as [x channel], but we send images back as [channel row col]?

[dchansen]

Presumably that would be entirely up to Gadgetron. Gadgetron Matlab is not based on ISMRMRD Matlab and would anyway need a compatibility layer for a new version of ISMRMRD.

[hansenms]

I think the sort answer is that an MRD version based on yardl would definitely break Gadgetron Matlab. Hopefully it would make it unnecessary in the long run if we support the intermediate data types too.

[kspaceKelvin]

Well, the existing situation is already confusing with MRD across languages, so it's a tricky starting point. The matrix_size field in the ImageHeader is in the order [x y z], which leads to the assumption that the data field itself would have the same shape (the addition of channels would be a mild surprise). However, as noted, the data field has shape data[channel, z, y, x] in C++ and Python and I think this has indeed been a source of confusion for some. In my Python example, I specifically note the dimensions and then permute it to [y, x, z, cha] before working on it.

It's worth pointing out that when creating an MRD image from Python using from_array, it assumes the data array will be [x y z cha] and will transpose the data for you:

>>> mrdImg = ismrmrd.Image.from_array(np.zeros((192,256,3,8)))
>>> mrdImg.data.shape
(8, 3, 256, 192)
>>> (mrdImg.getHead().matrix_size[0], mrdImg.getHead().matrix_size[1], mrdImg.getHead().matrix_size[2])
(192, 256, 3)

This does come with an associated deprecation warning:

/usr/local/lib/python3.10/site-packages/ismrmrd/image.py:191: PendingDeprecationWarning: The default behavior of this function is currently column-major which is inconsistent with numpy using row-major by default. In a future version this will be changed. Please switch to setting transpose in this function to false to switch to the new behavior.

Calling it with the transpose=False option would disable that behaviour:

>>> mrdImg = ismrmrd.Image.from_array(np.zeros((8,3,256,192)), transpose=False)
>>> mrdImg.data.shape
(8, 3, 256, 192)
>>> (mrdImg.getHead().matrix_size[0], mrdImg.getHead().matrix_size[1], mrdImg.getHead().matrix_size[2])
(192, 256, 3)

On the MATLAB side, the reference code in the ismrmrd repo doesn't actually have any code related to how the dimensions of image data is stored in +ismrmrd/Image.m. Whether you'd assume that implied [x y z cha] or the opposite probably depends on whether you started from the matrix_size header field or by looking at the serialization code.

In my fork, I mirrored the "transpose" of the Python from_array() call, as this also made sense to me when interpreting the matrix_size header field. This was prior to the Python deprecation warning being added too:

obj.head.matrix_size(1) = size(data,1);
obj.head.matrix_size(2) = size(data,2);
obj.head.matrix_size(3) = size(data,3);
obj.head.channels       = size(data,4);
obj.data                = data;

The Gadgetron fork the constructor for the Image class says self.data = data; (no transpose), and the deserialization code explicitly reshapes it into [channel z y x]. So there is already a discrepancy between the Gadgetron fork and my fork.

My personal opinion is that putting aside serialization and efficiency, accessing the data array as [x y z channel] makes the most intuitive sense -- the basic item is a 2D image, sometimes there's a 3rd dimension, and sometimes there's a multiple channels.

I agree with @ojosephs' point that [row col], i.e. [y x] as the first dimensions is the most convenient, especially when viewing the images. matplotlib's imshow also uses the [row col] ordering for displaying images. In my example code, I also explicitly convert to this shape:

# Reformat data to [y x z cha img], i.e. [row col] for the first two dimensions
data = data.transpose((3, 4, 2, 1, 0))

I also do this is MATLAB, though I guess I didn't put that one on GitHub yet. This is particularly important for integrating existing image processing code where it probably assumes [row col].

That said, I think the data inside the actual MRD class should be either [x y z channel] or [channel z y x], not [y x z channel]. My vote would be for [x y z channel], but I see the arguments for the opposite.

[KrisThielemans]

A lot of the comments here are related to ISMRMRD. We are using yardl for PETSIRD (and have no backwards compatibility issues). In any case, yardl could/will be used for many other projects. So I strongly suggest yardl-MATLAB handling should be logical in itself, and not be dictated by previous ISMRMRD conventions. (Of course, "logical" is virtually impossible with MATLAB in my opinion 😄.)

In MATLAB, swapping the first 2 dimensions is done in some functions, but not all. (I believe this is largely because MATLAB was designed for 2D, and then had to think how to extend it without breaking too many things.) An example is meshgrid (swapping first 2 coordinates) vs ndgrid (just using Fortran order).

IMO, the yardl-generated API should not do any swapping/permuting of data, for 2 reasons:

• MATLAB doesn't do any swapping either when reading/writing data. Indeed, if you read an nD-array in MATLAB and access it using linear indexing, the order is the same as when you do that in Python or C. (It is only one access it as an nD array, that you have to change the order of the indices).
• permutation is expensive, so shouldn't be done by the generated API, but by the user.

Thus, it is reasonable to establish a convention in yardl that Matlab N-dimensional arrays always have reversed dimensions with respect to the model.

I therefore 100% agree with this statement, and hence disagree with @ojosephs. In my opinion, using [y,x,z,col] would not be compatible with how MATLAB would read the data.

(The fact that imshow in MATLAB transposes is because their serialiser doesn't, and all 2D images therefore would have been displayed wrongly)

matrices

In yardl, there is no syntactic difference between ndarray and arrayOfMatrices in the example. Therefore, the generated API has to (de)serialise these objects in the same way. Given the above, it has to do that as an ndarray (without swapping of the first 2 dimensions) IMO. Similarly, there is no way in the example to know if the 2D example is intended as a matrix or as a 2D array. Therefore, it has to do the same as for an ndarray IMO.

If this is a problem for matrix calculations, then in my opinion the only way around that would be to create a yardl matrix type, which is handled differently than nD arrays. That seems easy enough to do, and clear as well. An arrayOfMatrices would then be of type matrix[count]. I'm not 100% sure this is necessary, but it is much better than doing different things based on a guess of what the intention of the user is.

vectors

This discussion hasn't touched on what to do for 1D arrays, which don't exist in MATLAB (at least, last time I checked). Is it a (1,n) or (n,1) array, which of course sometimes matters in MATLAB? My MATLAB is too rusty to remember what its own usual conventions are, but along the lines of a matrix type, yardl could have a columnvector and rowvector type.