API: decide on dimension ordering in spatial/kernel modules

[balbasty]

Here's the convention currently used in nitorch (copied from spatial.py):

---

NiTorch uses consistent ordering conventions throughout its API:

• We use abbreviations (B[atch], C[hannel], D[ephth], H[eight], W[idth]) to name dimensions.
• Tensors that represent series (resp. images or volumes) should always be ordered as (B, C, W) (resp. (B, C, H, W) or (B, C, D, H, W)).
• We use x, y, z to denote axis/coordinates along the (W, H, D) dimensions.
• Conversely, displacement or deformation fields are ordered as (B, H, D, W, C). There, the Channel dimension contains displacements or deformations along the x, y, z axes (i.e., W, H, D dimensions).
• Similarly, Jacobian fields are stored as (B, H, D, W, C, C). The second to last dimension corresponds to x, y, z components and the last dimension corresponds to derivatives along the x, y, z axes. I.e., jac[..., i, j] = d{u_i}/d{x_j}
• This means that we usually do not store deformations per imaging channel (they are assumed to lie in the same space).
• Arguments that relate to spatial dimensions are ordered as (x, y, z) or (W, H, D). This means: lattice dimensions, voxel sizes, affine transformation matrices, etc. It can be a little bit counterintuitive at first, as this order is opposite to the data storage order. However, it is consistant with the ordering of spatial components.

These conventions are mostly consistent with those used in PyTorch (conv, grid_sampler, etc.). However, some care must be taken when, e.g., convolving deformation fields or allocating tensors based on grid dimensions.

---

However, it is very counter-intuitive to sometime think in (D, H, W), and sometimes in (W, H, D). If we have to change this convention, it is better to do it very early in the project. After, it will be too late.

I think that what matters most is to be consistent with how nibabel loads volumes: in particular how rows/columns of the affine matrix map to dimensions of the loaded ndarray. We don't want users to have to permute dimensions of the loaded arrays.

@brudfors

[brudfors]

I agree this is an important one. At the moment I am permuting the nibabel array as .permute([2, 1, 0]) to make it compatible with nitorch.spatial. But if we don't permute the data, we would need to 'permute' the affine matrix right? So it seems to be an issue no matter what (in the case of many neuroimaing s/w and how they store data and affine transformations), or?

[balbasty]

For some reason, I thought that the dimension ordering from nibabel would be opposite to that of spm. Since they are the same, I agree that it is better to be consistent with nibabel.

I'll change everything so that all arguments are ordered as (D, H, W). This will break consistency with PyToch, but it is not a big deal. That will take a bit of time, I will create a branch to work on this.

[brudfors]

Sounds great! Will improve code readability too; e.g., in the superres the arrays holding x and y are permuted at the end of the algorithm, as to write to nifti with nibabel.