.. _coordmap-discussion:

########################################
Some discussion notes on coordinate maps
########################################

These notes contain some email discussion between Jonathan Taylor, Bertrand
Thirion and Gael Varoquaux about coordinate maps, coordinate systems and
transforms.

They are a little bit rough and undigested in their current form, but they might
be useful for background.

The code and discussion below mentions ideas like ``LPIImage``, ``XYZImage`` and
``AffineImage``.  These were image classes that constrained their coordinate
maps to have input and output axes in a particular order.  We eventually removed
these in favor of automated reordering of image axes on save, and explicit
reordering of images that needed known axis ordering.

.. some working notes

::

    import sympy
    i, j, k = sympy.symbols('i, j, k')
    np.dot(np.array([[0,0,1],[1,0,0],[0,1,0]]), np.array([i,j,k]))
    kij = CoordinateSystem('kij')
    ijk_to_kij = AffineTransform(ijk, kij, np.array([[0,0,1,0],[1,0,0,0],[0,1,0,0],[0,0,0,1]]))
    ijk_to_kij([i,j,k])
    kij = CoordinateSystem('kij')
    ijk_to_kij = AffineTransform(ijk, kij, np.array([[0,0,1,0],[1,0,0,0],[0,1,0,0],[0,0,0,1]]))
    ijk_to_kij([i,j,k])
    kij_to_RAS = compose(ijk_to_kij, ijk_to_RAS)
    kij_to_RAS = compose(ijk_to_RAS,ijk_to_kij)
    kij_to_RAS = compose(ijk_to_RAS,ijk_to_kij.inverse())
    kij_to_RAS
    kij = CoordinateSystem('kij')
    ijk_to_kij = AffineTransform(ijk, kij, np.array([[0,0,1,0],[1,0,0,0],[0,1,0,0],[0,0,0,1]]))
    # Check that it does the right permutation
    ijk_to_kij([i,j,k])
    # Yup, now let's try to make a kij_to_RAS transform
    # At first guess, we might try
    kij_to_RAS = compose(ijk_to_RAS,ijk_to_kij)
    # but we have a problem, we've asked for a composition that doesn't make sense
    kij_to_RAS = compose(ijk_to_RAS,ijk_to_kij.inverse())
    kij_to_RAS
    # check that things are working -- I should get the same value at i=20,j=30,k=40 for both mappings, only the arguments are reversed
    ijk_to_RAS([i,j,k])
    kij_to_RAS([k,i,j])
    another_kij_to_RAS = ijk_to_RAS.reordered_domain('kij')
    another_kij_to_RAS([k,i,j])
    # rather than finding the permutation matrix your self
    another_kij_to_RAS = ijk_to_RAS.reordered_domain('kij')
    another_kij_to_RAS([k,i,j])

    >>> ijk = CoordinateSystem('ijk', coord_dtype=np.array(sympy.Symbol('x')).dtype)
    >>> xyz = CoordinateSystem('xyz', coord_dtype=np.array(sympy.Symbol('x')).dtype)
    >>> x_start, y_start, z_start = [sympy.Symbol(s) for s in ['x_start', 'y_start', 'z_start']]
    >>> x_step, y_step, z_step = [sympy.Symbol(s) for s in ['x_step', 'y_step', 'z_step']]
    >>> i, j, k = [sympy.Symbol(s) for s in 'ijk']
    >>> T = np.array([[x_step,0,0,x_start],[0,y_step,0,y_start],[0,0,z_step,z_start],[0,0,0,1]])
    >>> T
    array([[x_step, 0, 0, x_start],
        [0, y_step, 0, y_start],
        [0, 0, z_step, z_start],
        [0, 0, 0, 1]], dtype=object)
    >>> A = AffineTransform(ijk, xyz, T)
    >>> A
    AffineTransform(
    function_domain=CoordinateSystem(coord_names=('i', 'j', 'k'), name='', coord_dtype=object),
    function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=object),
    affine=array([[x_step, 0, 0, x_start],
                    [0, y_step, 0, y_start],
                    [0, 0, z_step, z_start],
                    [0, 0, 0, 1]], dtype=object)
    )
    >>> A([i,j,k])
    array([x_start + i*x_step, y_start + j*y_step, z_start + k*z_step], dtype=object)
    >>> # this is another
    >>> A_kij = A.reordered_domain('kij')

    >>> A_kij
    AffineTransform(
    function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='', coord_dtype=object),
    function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=object),
    affine=array([[0, x_step, 0, x_start],
                    [0, 0, y_step, y_start],
                    [z_step, 0, 0, z_start],
                    [0.0, 0.0, 0.0, 1.0]], dtype=object)
    )
    >>>
    >>> A_kij([k,i,j])
    array([x_start + i*x_step, y_start + j*y_step, z_start + k*z_step], dtype=object)
                                                                                    >>> # let's look at another reordering
    >>> A_kij_yzx = A_kij.reordered_range('yzx')
    >>> A_kij_yzx
    AffineTransform(
    function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='', coord_dtype=object),
    function_range=CoordinateSystem(coord_names=('y', 'z', 'x'), name='', coord_dtype=object),
    affine=array([[0, 0, y_step, y_start],
                    [z_step, 0, 0, z_start],
                    [0, x_step, 0, x_start],
                    [0, 0, 0, 1.00000000000000]], dtype=object)
    )
    >>> A_kij_yzx([k,i,j])
    array([y_start + j*y_step, z_start + k*z_step, x_start + i*x_step], dtype=object)
    >>>

    class RASTransform(AffineTransform):
    """
    An AffineTransform with output, i.e. range:

    x: units of 1mm increasing from Right to Left
    y: units of 1mm increasing from Anterior to Posterior
    z:  units of 1mm increasing from Superior to Inferior
    """
    def reorder_range(self):
        raise ValueError('not allowed to reorder the "xyz" output coordinates')

    def to_LPS(self):
        from copy import copy
        return AffineTransform(copy(self.function_domain),
                                copy(self.function_range),
                                np.dot(np.diag([-1,-1,1,1], self.affine))

    class LPSTransform(AffineTransform):
    """
    An AffineTransform with output, i.e. range:

    x: units of 1mm increasing from Left to Right
    y: units of 1mm increasing from Posterior to Anterior
    z:  units of 1mm increasing from Inferior to Superior
    """
    def reorder_range(self):
        raise ValueError('not allowed to reorder the "xyz" output coordinates')


    def to_RAS(self):
        from copy import copy
        return AffineTransform(copy(self.function_domain),
                                copy(self.function_range),
                                np.dot(np.diag([-1,-1,1,1], self.affine)))

    class NeuroImage(Image):
    def __init__(self, data, affine, axis_names, world='world-RAS'):
        affine_transform = {'LPS':LPSTransform,
                            'RAS':RAITransform}[world])(axis_names[:3], "xyz", affine}
        ...

    LPIImage only forced it to be of one type.

Email #1
--------

Excuse the long email but I started writing, and then it started looking like documentation. I will put most of it into doc/users/coordinate_map.rst.


    Also, I am not sure what this means. The image is in LPI ordering, only
    if the reference frame of the world space it is pointing to is.


I am proposing we enforce the world space to have this frame of reference
to be explicit so that you could tell left from right on an image after calling xyz_ordered().


    If it is
    pointing to MNI152 (or Talairach), then x=Left to Right, y=Posterior to
    Anterior, and z=Inferior to Superior. If not, you are not in MNI152.
    Moreover, according to the FSL docs, the whole 'anatomical' versus
    'neurological' mess that I hear has been a long standing problem has
    nothing to do with the target frame of reference, but only with the way
    the data is stored.


I think the LPI designation simply specifies "x=Left to Right, y=Posterior to
Anterior, and z=Inferior to Superior" so any MNI152 or Tailarach would be in LPI
coordinates, that's all I'm trying to specify with the designation "LPI". If
MNI152 might imply a certain voxel size, then I would prefer not to use MNI152.

If there's a better colour for the bike shed, then I'll let someone else paint it, :)

This LPI specification actually makes a difference to the
"AffineImage/LPIImage.xyz_ordered" method. If, in the interest of being
explicit, we would enforce the direction of x,y,z in LPI/Neuro/AffineImage, then
the goal of having "xyz_ordered" return an image with an affine that has a
diagonal with positive entries, as in the AffineImage specification, means that
you might have to call

affine_image.get_data()[::-1,::-1] # or some other combination of flips

(i.e. you have to change how it is stored in memory).

The other way to return an diagonal affine with positive entries is to flip send
x to -x, y to -y, i.e. multiply the diagonal matrix by np.diag([-1,-1,1,1]) on
the left. But then your AffineImage would now have "x=Right to Left, y=Anterior
to Posterior" and we have lost the interpretation of x,y,z as LPI coordinates.

By being explicit about the direction of x,y,z we know that if the affine matrix
was diagonal and had a negative entry in the first position, then we know that
left and right were flipped when viewed with a command like::

    >>> pylab.imshow(image.get_data()[:,:,10])

Without specifying the direction of x,y,z we just don't know.

    You can of course create a new coordinate system describing, for instance
    the scanner space, where the first coordinate is not x, and the second
    not y, ... but I am not sure what this means: x, y, and z, as well as
    left or right, are just names. The only important information between two
    coordinate systems is the transform linking them.


The sentence:

"The only important information between two coordinate systems is the transform
linking them."

has, in one form or another, often been repeated in NiPy meetings, but no one
bothers to define the terms in this sentence.  So, I have to ask what is your
definition of "transform" and "coordinate system"?  I have a precise definition,
and the names are part of it.

Let's go through that sentence. Mathematically, if a transform is a function,
then a transform knows its domain and its range so it knows the what the
coordinate systems are. So yes, with transform defined as "function", if I give
you a transform between two coordinate systems (mathematical spaces of some
kind) the only important information about it is itself.

The problem is that, for a 4x4 matrix T, the python function

transform_function = lambda v: np.dot(T, np.hstack([v,1])[:3]

has a "duck-type" domain that knows nothing about image acquisition and a range inferred by numpy that knows nothing about LPI or MNI152.  The string "coord_sys" in AffineImage is meant to imply that its domain and range say it should be interpreted in some way, but it is not explicit in AffineImage.

(Somewhere around here, I start veering off into documentation.... sorry).

To me, a "coordinate system" is a basis for a vector space (sometimes you might
want transforms between integers but ignore them for now).  It's not even a
description of an affine subspace of a vector space, (see e.g.
http://en.wikipedia.org/wiki/Affine_transformation).  To describe such an affine
subspace, "coordinate system" would need one more piece of information, the
"constant" or "displacement" vector of the affine subspace.

Because it's a basis, each element in the basis can be identified by a name, so
the transform depends on the names because that's how I determine a "coordinate
system" and I need "coordinate systems" because they are what the domain and
range of my "transform" are going to be. For instance, this describes the range
"coordinate system" of a "transform" whose output is in LPI coordinates:

"x" = a unit vector of length 1mm pointing in the Left to Right direction
"y" = a unit vector of length 1mm pointing in the Posterior to Anterior direction
"z" = a unit vector of length 1mm pointing in the Inferior to Superior direction

OK, so that's my definition of "coordinate system" and the names are an
important part of it.

Now for the "transform" which I will restrict to be "affine transform". To me,
this is an affine function or transformation between two vector spaces (we're
not even considering affine transformations between affine spaces). I bring up
the distinction because generally affine transforms act on affine spaces rather
than vector spaces. A vector space is an affine subspace of itself with
"displacement" vector given by its origin, hence it is an affine space and so we
can define affine functions on vector spaces.

Because it is an affine function, the mathematical image of the domain under
this function is an affine subspace of its range (which is a vector space). The
"displacement" vector of this affine subspace is represented by the floats in b
where A,b = to_matvec(T) (once I have specified a basis for the range of this
function).

Since my "affine transform" is a function between two vector spaces, it should
have a domain that is a vector space, as well. For the "affine transform"
associated with an Image, this domain vector space has coordinates that can be
interpreted as array coordinates, or coordinates in a "data cube". Depending on
the acquisition parameters, these coordinates might have names like "phase",
"freq", "slice".

Now, I can encode all this information in a tuple: (T=a 4x4 matrix of floats
with bottom row [0,0,0,1], ('phase', 'freq', "slice"), ('x','y','z'))

>>> import numpy as np
>>> from nipy.core.api import CoordinateSystem, AffineTransform
>>> acquisition = ('phase', 'freq', 'slice')
>>> xyz_world = ('x','y','z')
>>> T = np.array([[2,0,0,-91.095],[0,2,0,-129.51],[0,0,2,-73.25],[0,0,0,1]])
>>> AffineTransform(CoordinateSystem(acquisition), CoordinateSystem(xyz_world), T)
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('phase', 'freq', 'slice'), name='', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=float64),
   affine=array([[   2.   ,    0.   ,    0.   ,  -91.095],
                 [   0.   ,    2.   ,    0.   , -129.51 ],
                 [   0.   ,    0.   ,    2.   ,  -73.25 ],
                 [   0.   ,    0.   ,    0.   ,    1.   ]])
)

The float64 appearing above is a way of specifying that the "coordinate systems"
are vector spaces over the real numbers, rather than, say the complex numbers.
It is specified as an optional argument to CoordinateSystem.

Compare this to the way a MINC file is described::

    jtaylo@ubuntu:~$ mincinfo data.mnc
    file: data.mnc
    image: signed__ short -32768 to 32767
    image dimensions: zspace yspace xspace
        dimension name         length         step        start
        --------------         ------         ----        -----
        zspace                     84            2       -73.25
        yspace                    114            2      -129.51
        xspace                     92            2      -91.095
    jtaylo@ubuntu:~$
    jtaylo@ubuntu:~$ mincheader data.mnc
    netcdf data {
    dimensions:
        zspace = 84 ;
        yspace = 114 ;
        xspace = 92 ;
    variables:
        double zspace ;
            zspace:varid = "MINC standard variable" ;
            zspace:vartype = "dimension____" ;
            zspace:version = "MINC Version    1.0" ;
            zspace:comments = "Z increases from patient inferior to superior" ;
            zspace:spacing = "regular__" ;
            zspace:alignment = "centre" ;
            zspace:step = 2. ;
            zspace:start = -73.25 ;
            zspace:units = "mm" ;
        double yspace ;
            yspace:varid = "MINC standard variable" ;
            yspace:vartype = "dimension____" ;
            yspace:version = "MINC Version    1.0" ;
            yspace:comments = "Y increases from patient posterior to anterior" ;
            yspace:spacing = "regular__" ;
            yspace:alignment = "centre" ;
            yspace:step = 2. ;
            yspace:start = -129.509994506836 ;
            yspace:units = "mm" ;
        double xspace ;
            xspace:varid = "MINC standard variable" ;
            xspace:vartype = "dimension____" ;
            xspace:version = "MINC Version    1.0" ;
            xspace:comments = "X increases from patient left to right" ;
            xspace:spacing = "regular__" ;
            xspace:alignment = "centre" ;
            xspace:step = 2. ;
            xspace:start = -91.0950012207031 ;
            xspace:units = "mm" ;
        short image(zspace, yspace, xspace) ;
            image:parent = "rootvariable" ;
            image:varid = "MINC standard variable" ;
            image:vartype = "group________" ;
            image:version = "MINC Version    1.0" ;
            image:complete = "true_" ;
            image:signtype = "signed__" ;
            image:valid_range = -32768., 32767. ;
            image:image-min = "--->image-min" ;
            image:image-max = "--->image-max" ;
        int rootvariable ;
            rootvariable:varid = "MINC standard variable" ;
            rootvariable:vartype = "group________" ;
            rootvariable:version = "MINC Version    1.0" ;
            rootvariable:parent = "" ;
            rootvariable:children = "image" ;
        double image-min ;
            image-min:varid = "MINC standard variable" ;
            image-min:vartype = "var_attribute" ;
            image-min:version = "MINC Version    1.0" ;
            image-min:_FillValue = 0. ;
            image-min:parent = "image" ;
        double image-max ;
            image-max:varid = "MINC standard variable" ;
            image-max:vartype = "var_attribute" ;
            image-max:version = "MINC Version    1.0" ;
            image-max:_FillValue = 1. ;
            image-max:parent = "image" ;
    data:

    zspace = 0 ;

    yspace = 0 ;

    xspace = 0 ;

    rootvariable = _ ;

    image-min = -50 ;

    image-max = 50 ;
    }

I like the MINC description, but the one thing missing in this file is the
ability to specify ('phase', 'freq', 'slice').  It may be possible to add it but
I'm not sure, it certainly can be added by adding a string to the header.  It
also mixes the definition of the basis with the affine transformation (look at
the output of mincheader which says that yspace has step 2). The NIFTI-1
standard allows limited possibilities to specify ('phase', 'freq', 'slice') this
with its dim_info byte but there are pulse sequences for which these names are
not appropriate.

One might ask: why bother making a "coordinate system" for the voxels. Well,
this is part of my definition of "affine transform".  More importantly, it
separates the notion of world axes ('x','y','z') and voxel indices
('i','j','k'). There is at least one use case, slice timing, a key step in the
fMRI pipeline, where we need to know which spatial axis is slice. One solution
would be to just add an attribute to AffineImage called "slice_axis" but then,
as Gael says, the possibilities for axis names are infinite, what if we want an
attribute for "group_axis"? AffineTransform provides an easy way to specify an
axis as "slice":

>>> unknown_acquisition = ('i','j','k')
>>> A = AffineTransform(CoordinateSystem(unknown_acquisition),
...                     CoordinateSystem(xyz_world), T)

After some deliberation, we find out that the third axis is slice...

>>> A.renamed_domain({'k':'slice'})
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j', 'slice'), name='', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=float64),
   affine=array([[   2.   ,    0.   ,    0.   ,  -91.095],
                 [   0.   ,    2.   ,    0.   , -129.51 ],
                 [   0.   ,    0.   ,    2.   ,  -73.25 ],
                 [   0.   ,    0.   ,    0.   ,    1.   ]])
)

Another question one might ask is: why bother allowing non-4x4 affine matrices
like:

>>> AffineTransform.from_params('ij', 'xyz', np.array([[2,3,1,0],[3,4,5,0],[7,9,3,1]]).T)
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j'), name='', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=float64),
   affine=array([[2., 3., 7.],
                 [3., 4., 9.],
                 [1., 5., 3.],
                 [0., 0., 1.]])
)

For one, it allows very clear specification of a 2-dimensional plane (i.e. a
2-dimensional affine subspace of some vector spce) called P, in, say, the  LPI
"coordinate system". Let's say we want the plane in LPI-world corresponding to
"j=30" for im above. (I guess that's coronal?)

Make an affine transform that maps (i,k) -> (i,30,k):

>>> j30 = AffineTransform(CoordinateSystem('ik'), CoordinateSystem('ijk'), np.array([[1,0,0],[0,0,30],[0,1,0],[0,0,1]]))
>>> j30
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'k'), name='', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('i', 'j', 'k'), name='', coord_dtype=float64),
   affine=array([[  1.,   0.,   0.],
                 [  0.,   0.,  30.],
                 [  0.,   1.,   0.],
                 [  0.,   0.,   1.]])
)

Its dtype is np.float since we didn't specify np.int in constructing the
CoordinateSystems:

>>> from nipy.core.api import compose
>>> j30_to_XYZ = compose(A, j30)
>>> j30_to_XYZ
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'k'), name='', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=float64),
   affine=array([[  2.   ,   0.   , -91.095],
                 [  0.   ,   0.   , -69.51 ],
                 [  0.   ,   2.   , -73.25 ],
                 [  0.   ,   0.   ,   1.   ]])
)

This could be used to resample any RAS Image on the coronal plane y=-69.51 with
voxels of size 2mm x 2mm starting at x=-91.095 and z=-73.25. Of course, this
doesn't seem like a very natural slice. The module
:mod:`nipy.core.reference.slices` has some convenience functions for specifying
slices.

>>> from nipy.core.reference.slices import yslice, bounding_box
>>> x_spec = ([-92,92], 93) # voxels of size 2 in x, starting at -92, ending at 92
>>> z_spec = ([-70,100], 86) # voxels of size 2 in z, starting at -70, ending at 100

When specifying a *y* slice - we have to know what "y" means.  In order for "y"
to have meaning, we need to specify the name of an output (range) space that has
a defined "y".  In this case we use MNI space:

>>> y70 = yslice(70, x_spec, z_spec, 'mni')
>>> y70
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i_x', 'i_z'), name='slice', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('mni-x=L->R', 'mni-y=P->A', 'mni-z=I->S'), name='mni', coord_dtype=float64),
   affine=array([[  2.,   0., -92.],
                 [  0.,   0.,  70.],
                 [  0.,   2., -70.],
                 [  0.,   0.,   1.]])
)

>>> x_lims, y_lims, z_lims = bounding_box(y70, (x_spec[1], z_spec[1]))
>>> assert np.all(x_lims == (-92, 92))
>>> assert np.all(y_lims == (70, 70))
>>> assert np.all(z_lims == (-70, 100))

Maybe these aren't things that "normal human beings" (to steal a quote from
Gael) can use, but they're explicit and they are tied to precise mathematical
objects.

Email #2
---------

I apologize again for the long emails, but I'm glad we. as a group, are having
this discussion electronically. Usually, our discussions of CoordinateMap begin
with Matthew standing in front of a white board with a marker and asking a
newcomer,

"Are you familiar with the notion of a transformation, say, from voxel to world?"

:)

Where they go after that really depends on the kind of day everyone's having...

:)

These last two emails also have the advantage that most of them can go right in
to doc/users/coordinate_map.rst.

    I agree with Gael that LPIImage is an obscure name.

OK. I already know that people often don't agree with names I choose, just ask
Matthew. :)

I just wanted to choose a name that is as explicit as possible. Since I'm
neither a neuroscientist nor an MRI physicist but a statistician, I have no idea
what it really means. I found it mentioned in this link below and John Ollinger
mentioned LPI in another email thread

http://afni.nimh.nih.gov/afni/community/board/read.php?f=1&i=9140&t=9140

I was suggesting we use a well-established term, apparently LPI is not
well-established. :)

Does LPS mean (left, posterior, superior)?  Doesn't that suggest that LPI means
(left, posterior, inferior) and RAI means (right, anterior, inferior)?  If so,
then good, now I know what LPI means and I'm not a neuroscientist or an MRI
physicist, :)

We can call the images RASImages, or at least let's call their AffineTransform
RASTransforms, or we could have NeuroImages that can only have RASTransforms or
LPSTransforms, NeuroTransform that have a property and NeuroImage raises an
exception like this::

    @property
    def world(self):
    return self.affine_transform.function_range

    if (self.world.name not in ['world-RAS', 'world-LPS'] or
        self.world.coord_names != ('x', 'y', 'z')):
        raise ValueError("the output space must be named one of "
                         "['world-RAS','world-LPS'] and "
                         "the axes must be ('x', 'y', 'z')")

    _doc['world'] = "World space, one of ['world-RAS', 'world-LPS']. If it is 'world-LPS', then x increases from patient's left to right, y increases posterior to anterior, z increases superior to inferior. If it is 'world-RAS' then x increases patient's right to left, y increases posterior to anterior, z increases superior to inferior."

I completely advocate any responsibility for deciding which acronym to choose,
someone who can use rope can just change every lpi/LPI to ras/RAS I just want it
explicit.  I also want some version of these phrases "x increases from patient's
right to left", "y increases from posterior to anterior", "z increases from
superior to inferior" somewhere in a docstring for RAS/LPSTransform (see why I
feel that "increasing vs. decreasing" is important below).

I want the name and its docstring to scream at you what it represents so there
is no discussion like on the AFNI list where users are not sure which output of
which program (in AFNI) should be flipped (see the other emails in the thread).
It should be a subclass of AffineTransform because it has restrictions: namely,
its range is 'xyz'  and "xy" can be interpreted in of two ways either RAS or
LPS). You can represent any other version of RAS/LPS or (whatever colour your
bike shed is, :)) with the same class, it just may have negative values on the
diagonal. If it has some rotation applied, then it becomes pretty hard (at least
for me) to decide if it's RAS or LPS from the 4x4 matrix of floats. I can't even
tell you now when I look at the FIAC data which way left and right go unless I
ask Matthew.

    For background, you may want to look at what Gordon Kindlmann did for
    nrrd format where you can declare the space in which your orientation
    information and other transforms should be interpreted:

    http://teem.sourceforge.net/nrrd/format.html#space

    Or, if that's too flexible for you, you could adopt a standard space.

    ITK chose LPS to match DICOM.

    For slicer, like nifti, we chose RAS

It may be that there is well-established convention for this, but then why does
ITK say DICOM=LPS and AFNI say DICOM=RAI?  At least MINC is explicit. I favor
making it as precise as MINC does.

That AFNI discussion I pointed to uses the pairing RAI/DICOM and LPI/SPM.  This
discrepancy suggests there's some disagreement between using the letters to name
the system and whether they mean increasing or decreasing. My guess is that
LPI=RAS based on ITK/AFNI's identifications of LPS=DICOM=RAI. But I can't tell
if the acronym LPI means "x is increasing L to R, y increasing from P to A, z in
increasing from I to S" which would be equivalent to RAS meaning "x decreasing
from R to L, y decreasing from A to P, z is decreasing from S to I". That is, I
can't tell from the acronyms which of LPI or RAS is using "increasing" and which
is "decreasing", i.e. they could have flipped everything so that LPI means "x is
decreasing L to R, y is decreasing P to A, z is decreasing I to S" and RAS means
"x is increasing R to L, y is increasing A to P, z is increasing S to I".

To add more confusion to the mix, the acronym doesn't say if it is the patient's
left to right or the technician looking at him, :) For this, I'm sure there's a
standard answer, and it's likely the patient, but heck, I'm just a statistician
so I don't know the answer.

    (every volume has an ijkToRAS affine transform).  We convert to/from LPS
    when calling ITK code, e.g., for I/O.

How much clearer can you express "ijkToRAS" or "convert to/from LPS" than
something like this:

>>> T = np.array([[2,0,0,-91.095],[0,2,0,-129.51],[0,0,2,-73.25],[0,0,0,1]])
>>> ijk = CoordinateSystem('ijk', 'voxel')
>>> RAS = CoordinateSystem('xyz', 'world-RAS')
>>> ijk_to_RAS = AffineTransform(ijk, RAS, T)
>>> ijk_to_RAS
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j', 'k'), name='voxel', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='world-RAS', coord_dtype=float64),
   affine=array([[   2.   ,    0.   ,    0.   ,  -91.095],
                 [   0.   ,    2.   ,    0.   , -129.51 ],
                 [   0.   ,    0.   ,    2.   ,  -73.25 ],
                 [   0.   ,    0.   ,    0.   ,    1.   ]])
)

>>> LPS = CoordinateSystem('xyz', 'world-LPS')
>>> RAS_to_LPS = AffineTransform(RAS, LPS, np.diag([-1,-1,1,1]))
>>> ijk_to_LPS = compose(RAS_to_LPS, ijk_to_RAS)
>>> RAS_to_LPS
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('x', 'y', 'z'), name='world-RAS', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='world-LPS', coord_dtype=float64),
   affine=array([[-1.,  0.,  0.,  0.],
                 [ 0., -1.,  0.,  0.],
                 [ 0.,  0.,  1.,  0.],
                 [ 0.,  0.,  0.,  1.]])
)
>>> ijk_to_LPS
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j', 'k'), name='voxel', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='world-LPS', coord_dtype=float64),
   affine=array([[  -2.   ,    0.   ,    0.   ,   91.095],
                 [   0.   ,   -2.   ,    0.   ,  129.51 ],
                 [   0.   ,    0.   ,    2.   ,  -73.25 ],
                 [   0.   ,    0.   ,    0.   ,    1.   ]])
)

Of course, we shouldn't rely on the names ijk_to_RAS to know that it is an
ijk_to_RAS transform, that's why they're in the AffineTransform. I don't think
any one wants an attribute named "ijk_to_RAS" for AffineImage/Image/LPIImage.

The other problem that LPI/RAI/AffineTransform addresses is that someday you
might want to transpose the data in your array and still have what you would
call an "image". AffineImage allows this explicitly because there is no
identifier for the domain of the AffineTransform (the attribute name "coord_sys"
implies that it refers to either the domain or the range but not both). (Even
those who share the sentiment that "everything that is important about the
linking between two coordinate systems is contained in the transform"
acknowledge there are two coordinate systems :))

Once you've transposed the array, say

>>> data = np.random.normal(size=(10, 12, 14)) # original array
>>> newdata = data.transpose([2,0,1])

You shouldn't use something called "ijk_to_RAS" or "ijk_to_LPS" transform.
Rather, you should use a "kij_to_RAS" or "kij_to_LPS" transform.

>>> ijk = CoordinateSystem('ijk', 'voxel')
>>> kij = CoordinateSystem('kij', 'voxel')
>>> ijk_to_kij = AffineTransform(ijk, kij, np.array([[0,0,1,0],[1,0,0,0],[0,1,0,0],[0,0,0,1]]))

Check that it does the right permutation

>>> i, j, k = 10., 20., 40
>>> ijk_to_kij([i, j, k])
array([40., 10., 20.])

Yup, now let's try to make a kij_to_RAS transform

At first guess, we might try

>>> kij_to_RAS = compose(ijk_to_RAS, ijk_to_kij)
Traceback (most recent call last):
    ...
ValueError: domains and ranges don't match up correctly

We have a problem, we've asked for a composition that doesn't make sense.

If you're good with permutation matrices, you wouldn't have to call "compose"
above and you can just do matrix multiplication.  But here the name of the
function tells you that yes, you should do the inverse: "ijk_to_kij" says that
the range are "kij" values, but to get a "transform" for your data in "kij" it
should have a domain that is "kij".

The call to compose raised an exception because it saw you were trying to
compose a function with domain="ijk" and range="kji" with a function (on its
left) having domain="ijk" and range "kji". This composition just doesn't make
sense so it raises an exception.

>>> kij_to_ijk = ijk_to_kij.inverse()
>>> kij_to_RAS = compose(ijk_to_RAS, kij_to_ijk)
>>> kij_to_RAS
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='voxel', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='world-RAS', coord_dtype=float64),
   affine=array([[   0.   ,    2.   ,    0.   ,  -91.095],
                 [   0.   ,    0.   ,    2.   , -129.51 ],
                 [   2.   ,    0.   ,    0.   ,  -73.25 ],
                 [   0.   ,    0.   ,    0.   ,    1.   ]])
)


>>> ijk_to_RAS([i,j,k])
array([-71.095, -89.51 ,   6.75 ])
>>> kij_to_RAS([k,i,j])
array([-71.095, -89.51 ,   6.75 ])

We also shouldn't have to rely on the names of the AffineTransforms, i.e.
ijk_to_RAS,  to remember what's what (in typing this example, I mixed up kij and
kji many times). The objects ijk_to_RAS, kij_to_RAS represent the same "affine
transform", as evidenced by their output above.  There are lots of
representations of the same "affine transform": (6=permutations of
i,j,k)*(6=permutations of x,y,z)=36 matrices for one "affine transform".

If we throw in ambiguity about the sign in front of the output, there are
36*(8=2^3 possible flips of the x,y,z)=288 matrices possible but there are only
really 8 different "affine transforms". If you force the order of the range to
be "xyz" then there are 6*8=48 different matrices possible, again only
specifying 8 different "affine transforms". For AffineImage, if we were to allow
both "LPS" and "RAS" this means two flips are allowed, namely either
"LPS"=[-1,-1,1] or "RAS"=[1,1,1], so there are 6*2=12 possible matrices to
represent 2 different "affine transforms".

Here's another example that uses sympy to show what's going on in the 4x4 matrix
as you reorder the 'ijk' and the 'RAS'. (Note that this code won't work in
general because I had temporarily disabled a check in CoordinateSystem that
enforced the dtype of the array to be a builtin scalar dtype for sanity's sake).
To me, each of A, A_kij and A_kij_yzx below represent the same "transform"
because if I substitute i=30, j=40, k=50 and I know the order of the 'xyz' in the
output then they will all give me the same answer.

>>> import sympy
>>> ijk = CoordinateSystem('ijk', coord_dtype=np.array(sympy.Symbol('x')).dtype)
>>> xyz = CoordinateSystem('xyz', coord_dtype=np.array(sympy.Symbol('x')).dtype)
>>> x_start, y_start, z_start = [sympy.Symbol(s) for s in ['x_start', 'y_start', 'z_start']]
>>> x_step, y_step, z_step = [sympy.Symbol(s) for s in ['x_step', 'y_step', 'z_step']]
>>> i, j, k = [sympy.Symbol(s) for s in 'ijk']
>>> T = np.array([[x_step,0,0,x_start],[0,y_step,0,y_start],[0,0,z_step,z_start],[0,0,0,1]])
>>> T
array([[x_step, 0, 0, x_start],
       [0, y_step, 0, y_start],
       [0, 0, z_step, z_start],
       [0, 0, 0, 1]], dtype=object)
>>> A = AffineTransform(ijk, xyz, T)
>>> A
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j', 'k'), name='', coord_dtype=object),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=object),
   affine=array([[x_step, 0, 0, x_start],
                 [0, y_step, 0, y_start],
                 [0, 0, z_step, z_start],
                 [0, 0, 0, 1]], dtype=object)
)
>>> A([i,j,k]) == [x_start + i*x_step, y_start + j*y_step, z_start + k*z_step]
array([ True,  True,  True])

This is another

>>> A_kij = A.reordered_domain('kij')
>>> A_kij
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='', coord_dtype=object),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=object),
   affine=array([[0, 1.0*x_step, 0, 1.0*x_start],
                 [0, 0, 1.0*y_step, 1.0*y_start],
                 [1.0*z_step, 0, 0, 1.0*z_start],
                 [0.0, 0.0, 0.0, 1.0]], dtype=object)
)
>>> A_kij([k,i,j])
array([1.0*i*x_step + 1.0*x_start, 1.0*j*y_step + 1.0*y_start,
       1.0*k*z_step + 1.0*z_start], dtype=object)

Let's look at another reordering:

>>> A_kij_yzx = A_kij.reordered_range('yzx')
>>> A_kij_yzx
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='', coord_dtype=object),
   function_range=CoordinateSystem(coord_names=('y', 'z', 'x'), name='', coord_dtype=object),
   affine=array([[0, 0, 1.0*y_step, 1.0*y_start],
                 [1.0*z_step, 0, 0, 1.0*z_start],
                 [0, 1.0*x_step, 0, 1.0*x_start],
                 [0, 0, 0, 1.00000000000000]], dtype=object)
)
>>> A_kij_yzx([k,i,j])
array([1.0*j*y_step + 1.0*y_start, 1.0*k*z_step + 1.0*z_start,
       1.0*i*x_step + 1.0*x_start], dtype=object)

>>> A_kij
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('k', 'i', 'j'), name='', coord_dtype=object),
   function_range=CoordinateSystem(coord_names=('x', 'y', 'z'), name='', coord_dtype=object),
   affine=array([[0, 1.0*x_step, 0, 1.0*x_start],
                 [0, 0, 1.0*y_step, 1.0*y_start],
                 [1.0*z_step, 0, 0, 1.0*z_start],
                 [0.0, 0.0, 0.0, 1.0]], dtype=object)
)

>>> from nipy.core.reference.coordinate_map import equivalent
>>> equivalent(A_kij, A)
True
>>> equivalent(A_kij, A_kij_yzx)
True
