Understanding voxel and real world mappings¶
Voxel coordinates and real-world coordinates¶
A point can be represented by coordinates relative to specified axes. coordinates are (almost always) numbers - see coordinate systems
For example, a map grid reference gives a coordinate (a pair of
numbers) to a point on the map. The numbers give the respective
positions on the horizontal (x) and vertical (y) axes of the
map.
A coordinate system is defined by a set of axes. In the example
above, the axes are the x and y axes. Axes for coordinates
are usually orthogonal - for example, moving one unit up on the x
axis on the map causes no change in the y coordinate - because
the axes are at 90 degrees.
In this discussion we’ll concentrate on the three dimensional case. Having three dimensions means that we have a three axis coordinate system, and coordinates have three values. The meaning of the values depend on what the axes are.
Voxel coordinates¶
Array indexing is one example of using a coordinate system. Let’s say we have a three dimensional array:
A = np.arange(24).reshape((2,3,4))
The value 0 is at array coordinate 0,0,0:
assert A[0,0,0] == 0
and the value 23 is at array coordinate 1,2,3:
assert A[1,2,3] == 23
(remembering python’s zero-based indexing). If we now say that our array is a 3D volume element array - an array of voxels, then the array coordinate is also a voxel coordinate.
If we want to use numpy to index our array, then we need integer
voxel coordinates, but if we use a resampling scheme, we can also
imagine non-integer voxel coordinates for A, such as
(0.6,1.2,1.9), and we could use resampling to estimate the value
at such a coordinate, given the actual data in the surrounding
(integer) points.
Array / voxel coordinates refer to the array axes. Without any
further information, they do not tell us about where the point is in
the real world - the world we can measure with a ruler. We refer to
array / voxel coordinates with indices i, j, k, where i is the
first value in the 3 value coordinate tuple. For example, if array /
voxel point (1,2,3) has i=1, j=2, k=3. We’ll be careful only
to use i, j, k rather than x, y, z, because we are going to
use x, y, z to refer to real-world coordinates.
Real-world coordinates¶
Real-world coordinates are coordinates where the values refer to real-world axes. A real-world axis is an axis that refers to some real physical space, like low to high position in an MRI scanner, or the position in terms of the subject’s head.
Here we’ll use the usual neuroimaging convention, and that is to label our axes relative to the subject’s head:
xhas negative values for left and positive values for right
yhas negative values for posterior (back of head) and positive values for anterior (front of head)
zhas negative values for the inferior (towards the neck) and positive values for superior (towards the highest point of the head, when standing)
Image index ordering¶
Background¶
In general, images - and in particular NIfTI format images, are ordered in memory with the X dimension changing fastest, and the Z dimension changing slowest.
Numpy has two different ways of indexing arrays in memory, C and fortran. With C index ordering, the first index into an array indexes the slowest changing dimension, and the last indexes the fastest changing dimension. With fortran ordering, the first index refers to the fastest changing dimension - X in the case of the image mentioned above.
C is the default index ordering for arrays in Numpy.
For example, let’s imagine that we have a binary block of 3D image data, in standard NIfTI / Analyze format, with the X dimension changing fastest, called my.img, containing Float32 data. Then we memory map it:
img_arr = memmap('my.img', dtype=float32)
When we index this new array, the first index indexes the Z dimension, and the third indexes X. For example, if I want a voxel X=3, Y=10, Z=20 (zero-based), I have to get this from the array with:
img_arr[20, 10, 3]
The problem¶
Most potential users of NiPy are likely to have experience of using image arrays in Matlab and SPM. Matlab uses Fortran index ordering. For fortran, the first index is the fastest changing, and the last is the slowest-changing. For example, here is how to get voxel X=3, Y=10, Z=20 (zero-based) using SPM in Matlab:
img_arr = spm_read_vols(spm_vol('my.img'));
img_arr(4, 11, 21) % matlab indexing is one-based
This ordering fits better with the way that we talk about coordinates in functional imaging, as we invariably use XYZ ordered coordinates in papers. It is possible to do the same in numpy, by specifying that the image should have fortran index ordering:
img_arr = memmap('my.img', dtype=float32, order='F')
img_arr[3, 10, 20]
The proposal¶
Change the default ordering of image arrays to fortran, in order to allow XYZ index ordering. So, change the access to the image array in the image class so that, to get the voxel at X=3, Y=10, Z=20 (zero-based):
img = Image('my.img')
img[3, 10, 20]
instead of the current situation, which requires:
img = Image('my.img')
img[20, 10, 3]
Summary of discussion¶
For:
Fortran index ordering is more intuitive for functional imaging because of conventional XYZ ordering of spatial coordinates, and Fortran index ordering in packages such as Matlab
Indexing into a raw array is fast, and common in lower-level applications, so it would be useful to implement the more intuitive XYZ ordering at this level rather than via interpolators (see below)
Standardizing to one index ordering (XYZ) would mean users would not have to think about the arrangement of the image in memory
Against:
C index ordering is more familiar to C users
C index ordering is the default in numpy
XYZ ordering can be implemented by wrapping by an interpolator
Potential problems¶
Performance penalties¶
KY commented:
This seems like a good idea to me but I have no knowledge of numpy
internals (and even less than none after the numeric/numarray
integration). Does anyone know if this will (or definitely will not)
incur any kind of obvious performance penalties re. array operations
(sans arcane problems like stride issues in huge arrays)?
MB replied:
Note that, we are not proposing to change the memory layout of the image, which is fixed by the image format in e.g NIfTI, but only to index it XYZ instead of ZYX. As far as I am aware, there are no significant performance differences between:
img_arr = memmap('my.img', dtype=float32, order='C') img_arr[5,4,3]and:
img_arr = memmap('my.img', dtype=float32, order='F') img_arr[3,4,5]Happy to be corrected though.
Clash between default ordering of numpy arrays and nipy images¶
C index ordering is the default in numpy, and using fortran ordering for images might be confusing in some circumstances. Consider for example:
img_obj = Image(‘my.img’) # Where the Image class has been changed to implement Fortran ordering first_z_slice = img_obj[…,0] # returns a Z slice
img_arr = memmap(‘my.img’, dtype=float32) # C ordering, the numpy default img_obj = Image(img_arr) first_z_slice = img_obj[…,0] # in fact returns an X slice
I suppose that we could check that arrays are fortran index ordered in the Image __init__ routine.
An alternative proposal - XYZ ordering of output coordinates¶
JT: Another thought, that is a compromise between the XYZ coordinates and Fortran ordering.
To me, having worked mostly with C-type arrays, when I index an array I think in C terms. But, the Image objects have the “warp” attached to them, which describes the output coordinates. We could insist that the output coordinates are XYZT (or make this an option). So, for instance, if the 4x4 transform was the identity, the following two calls would give something like:
interp = interpolator(img)
img[3,4,5] == interp(5,4,3)
This way, users would be sure in the interpolator of the order of the coordinates, but users who want access to the array would know that they would be using the array order on disk…
I see that a lot of users will want to think of the first coordinate as “x”, but depending on the sampling the [0] slice of img may be the leftmost or the rightmost. To find out which is which, users will have to look at the 4x4 transform (or equivalently the start and the step). So just knowing the first array coordinate is the “x” coordinate still misses some information, all of which is contained in the transform.
MB replied:
I agree that the output coordinates are very important - and I think we all agree that this should be XYZ(T)?
For the raw array indices - it is very common for people to want to do things to the raw image array - the quickstart examples containing a few - and you usually don’t care about which end of X is left in that situation, only which spatial etc dimension the index refers to.