Python geometry order of axes

[genusn]

This issue is related to #157, #189

I would like to confirm the order of x, y, z for geometry. In othrer words, I would like to know which is expected, for example, nVoxel=np.array([nVoxelZ, nVoxelY, nVoxelX]) or nVoxel=np.array([nVoxelX, nVoxelY, nVoxelZ])?

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import tigre
from tigre.demos.Test_data import data_loader

nangles=30
nVoxelX=256
nVoxelY=128
nVoxelZ=64
print("nVoxelX, nVoxelY, nVoxelZ = {}, {}, {}".format(nVoxelX, nVoxelY, nVoxelZ))
geo=tigre.geometry(mode='cone', nVoxel=np.array([nVoxelZ, nVoxelY, nVoxelX]), default=True)
print("geo.nVoxel = {}".format(geo.nVoxel))
angles=np.linspace(0,2*np.pi-2*np.pi/nangles,nangles)
head=data_loader.load_head_phantom(geo.nVoxel)

listIdxZ = [0, nVoxelZ*3//8, nVoxelZ//2, nVoxelZ*5//8, nVoxelZ-1]

fig, axes = plt.subplots(len(listIdxZ),1, figsize=(5, 8/0.97), sharey=True)
fig.suptitle('Inisotropic volume, Head, geo.nVoxel = {}'.format(geo.nVoxel))

for idx in range(len(listIdxZ)):
    idxZ = listIdxZ[idx]
    axis = axes[idx]
    axis.set_ylabel("z={}".format(idxZ))
    mappable = axis.imshow(head[idxZ])
    divider = make_axes_locatable(axis)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    fig.colorbar(mappable, cax=cax)

plt.tight_layout(rect=[0,0,1,0.97])
fig.savefig("InosotropicVolumeHead_Python.png")

results

and if modified to nVoxel=np.array([nVoxelX, nVoxelY, nVoxelZ] is passed to tigre.geometry, the result is

This means that data_loader.load_head_phantom() assumes geo.nVoxel is ZYX order. Is this correct?

[AnderBiguri]

This means that data_loader.load_head_phantom() assumes geo.nVoxel is ZYX order. Is this correct?

Yes I think so.

To be honest, this is something that was put together without too much though and never properly finished, so how we want the python version to be is something that we can discuss now, ignoring what TIGRE currently has.

My guts tell me that the arrays need to be ZYX (or ZXY?, the other one seems better), but that when we declare it in the geometry, we need to do X,Y,Z, because it feels easier. But maybe that is my matlab talking, and we should change the geo class to get the shapes and sizes in the same order as the array indices?

What do you think would be best?

[genusn]

the arrays need to be ZYX

I agree with you. To make the things easier (for user and developer), I think that x should be the fastest and z be the slowest variable on the linear memory space. If we accept it, we need also to accept img[z,y,x] for accessing the point at (x,y,z).

we declare it in the geometry, we need to do X,Y,Z, because it feels easier

Half agree, half disagree... very difficult decision, however, if we accept img[z,y,x], I think zyx would be better because all users have to remember becomes just that the order of components in python is z, y, x.

we should change the geo class to get the shapes and sizes in the same order as the array indices?

If we do so, I think that necessary is to change _types.pxd. The Geometry class is only a collection of arrays and no item is assigned to any of X,Y,Z, whereas the correspondence to X, Y, Z are defined in convert_to_c_geometry() of _types.pxd. I doubt Geometry.convert_contig_mode() is necessary because it is not the python geo object itself but C version copied (and converted) from it that is passed to the C world.

[AnderBiguri]

Half agree, half disagree... very difficult decision, however, if we accept img[z,y,x], I think zyx would be better because all users have to remember becomes just that the order of components in python is z, y, x.

But if we decide to go zyx for the geo arrays, then we should very strongly document this, as while it fits with the array shape, it does not fit with our human brains! An alternative can also be to have named variables, such as nVoxelZ etc, instead of arrays, maybe have Geometry class with getter/setters? Just brainstorming.

I doubt Geometry.convert_contig_mode() is necessary

plausible, I had less say than I should have (because I delegated, my fault!) into the python code, so its possible that its not needed. I don't remember why its there.... Basically, if you see anything that you feel should not be there, you are possibly correct :) feel free to change it.

[genusn]

Importance of documentation. 100% agree. Introducing getter/setters is a good idea especially for roll-pitch-yaw!

I found pages about multi-dimensional images in numpy. I think they are encouraging us to choose zyx.

1. https://numpy.org/doc/stable/reference/internals.html#multidimensional-array-indexing-order-issues
2. https://scikit-image.org/docs/dev/user_guide/numpy_images.html#coordinate-conventions

[AnderBiguri]

Ah, reading the second link, I see all the confusions that I am having.

It says that the convention is:

3D grayscale | (pln, row, col)

and python being a language with row-major ordering, this means that if unrolled to 1D, its the second axis that will appear first, then the third axis, finally the first axis.

MATLAB (and the CUDA code in TIGRE) is designed so that the 3D grayscale image is (col, row, pln), because MATLAB being column-major, the first axis is the one unrolled, then second, then third. So, for both memory data to be the same, a matrix in MATLAB would be XYZ, then the one in python would be ZXY. This would mean that the CUDA code needs changing.

Now, when this code was being developed, I was under the impression that we were using the same memory layout for MATLAB and python, i.e. python was ZXY. But if we are now doing ZYX, then all the CUDA problems are fixed, as by changing the axis-to-memory map, we "accidentally" accommodate for code that was intended for column-major arrays.

I am happy with this, ZYX is good. However, I have a worry: python does indeed expect that the second axis is the fastest one, coalescent memory. This means that it will write code assuming that, and when you break the assumption, all numpy operations can actually go quite slower (e.g. try multiplying two matrices, then transposing one and multiplying them again, numpy transpose does not change memory layout and thus makes operations quite slower). My worry is that considering all this is not how I was expecting, the code in algorithms etc may not be optimized for ZYX but for ZXY instead. I have no idea, we'll need to do some profiling at some point, I guess.

[genusn]

3D grayscale | (pln, row, col) python does indeed expect that the second axis is the fastest one, coalescent memory

I think that python expects that pln is the slowest and col is the fastest, and row (the second axis) is the second fastest/slowest.

As long as I have tested the current TIGRE, the comparison of MATLAB and Python is as follows:

3D image

As for the point of (X,Y,Z)

	MATLAB	Python
Voxel value	mat(iX,iY,iZ)	mat[iZ,iY,iX]
Axial slice	mat(:,:,iZ)	mat[iZ]
Memory layout	iX+nVoxelX * (iY + nVoxelY * iZ )	Same as left

Projection image

As for the point of (U,V) of T-th Projection

	MATLAB	Python
Pixel value	~~mat(iU,iV,iT)~~ mat(iV,iU,iT)	mat[iT,iV,iU]
Projection	mat(:,:,iT)	mat[iT]
Sinogram	~~mat(:,iV,:)~~ mat(iV,:,:)	mat[:,iV,:]
Memory layout	iV+nDetectorV * (iU + nDetectorU * iT )	iU+nDetectorU * (iV + nDetectorV * iT )

[AnderBiguri]

I need to do more tests.... but I am still unconvinced, MATLAB and python do, in fact, store memory in a different order, so if the CUDA code works, there must be somewhere we are changing the convention. I can't see where, if what you say is correct (which you prove correct)

This can be tested simply with:

import numpy as np

a=np.array([0,1,2,3,4,5,6,7])
b=a.reshape([2,2,2])

print(b.ravel(order = 'C')) # Python memory layout
print(b.ravel(order = 'F')) # MATLAB memory layout

As when we get to CUDA, all the memory layout is what we get, without a copy, for it to work in both, we must have the python one already by convention in the MATLAB layout. Let me think a bit, this seems more of my brain issue than code issue.

By the way, MATLAB uses mat(iV,iU,iT) for projections. Not sure if you where just randomly naming the indices, but the first index in TIGRE projections is the one corresponding to V axes. Not ideal, but too late to change.

[AnderBiguri]

OK, sorry, I understood it!

MATLAB indexing goes <row,col,plt>, and python goes <plt,row,col>. The swap of memory layout happening because row/column major gets compensated because in MATLAB we call "X" the first axis, i.e. the one corresponding to the row of the matrix while on python, its instead the col. That swap offsets the one in the memory layout.

Sorry, my brain needed to process this :) you are correct, of course.

[genusn]

By the way, MATLAB uses mat(iV,iU,iT) for projections.

I am sorry. You are right. Corrected the table above.

[AnderBiguri]

@genusn no worries, its a reasonable mistake, it should be as you wrote logically, bur for some reason I decided to not do that back then. I think thats a mistake from my side, but likely too big of a task to change now.

[genusn]

By the way, MATLAB uses mat(iV,iU,iT) for projections.

I am sorry. You are right. Corrected the table above.

@AnderBiguri ... so in MATLAB, despite the pixel value is mat(iV,iU,iT), nDetector has to be declared in the U-V (not V-U) order. One can check this by running the following code:

clear;
close all;
%% Geometry
nVoxelX=256;
nVoxelY=128;
nVoxelZ=64;
nDetectorU=250;
nDetectorV=125;
geo=defaultGeometry('nVoxel',[nVoxelX;nVoxelY;nVoxelZ]);
geo.nDetector=[nDetectorU; nDetectorV];	%  U-V order
geo.dDetector=[0.8*2; 0.8*2]; 
geo.sDetector=geo.nDetector.*geo.dDetector;

angles=linspace(0,2*pi,100);
head=headPhantom(geo.nVoxel);
projections=Ax(head,geo,angles,'interpolated');

projections(int32(nDetectorV*2/4-5 : nDetectorV) ...
   , int32(nDetectorU*1/4-5 : nDetectorU*1/4+5) ...
    , :) = 30;     % V-U-T order.

vol=Atb(projections, geo,angles);

plotProj(projections, angles)

Because of this, MLEM stops for non-square detectors. Modifying the shape of the weight matrix in https://github.com/CERN/TIGRE/blob/d707f0c866b8a02c7038d397408e5d4213b9e8cd/MATLAB/Algorithms/MLEM.m#L26 to [geo.nDetector(2) geo.nDetector(1) numel(angles)] was a solution to make MLEM.m work.

[AnderBiguri]

@genusn ah that is correct. I'll try to modify that. Admittedly, MLEM is a bit out of place and never cared too much about it.... CT data can not really be modeled very well as a Poisson distribution, and even in applications where the data can be (such as Positron Emission Tomography) MLEM is an obsolete algorithm because it needs hundreds of iterations to converge to the ML solution... So I always considered removing it, but at the same time, I think that its better its there and people can try it.

Just excusing the little care I gave MLEM :) I'll get that fixed.

[genusn]

Revised the table above. The xVoxel/xDetector rows are added.

3D image

As for the point of (X,Y,Z)

	MATLAB	Python
geo.xVoxel	[X;Y;Z]	[Z,Y,X]
Voxel value	mat(iX,iY,iZ)	mat[iZ,iY,iX]
Axial slice	mat(:,:,iZ)	mat[iZ]
Memory layout	iX+nVoxelX * (iY + nVoxelY * iZ )	Same as left

Projection image

As for the point of (U,V) of T-th Projection

	MATLAB	Python
geo.xDetector	[U;V]	[V,U]
Pixel value	mat(iV,iU,iT)	mat[iT,iV,iU]
Projection	mat(:,:,iT)	mat[iT]
Sinogram	mat(iV,:,:)	mat[:,iV,:]
Memory layout	iV+nDetectorV * (iU + nDetectorU * iT )	iU+nDetectorU * (iV + nDetectorV * iT )