RF,ENH: Refactor Nifti class - it now pre-calculates affines between voxel,

fsl, and world coordinate systems. New getAffine method allows you to get an
affine between any of the coordinate systems

fsl/data/image.py

@@ -135,13 +135,32 @@ class Nifti(notifier.Notifier, meta.Meta):
     methods.
 
 
-    **The affine transformation**
+    **Affine transformations**
 
 
-    The :meth:`voxToWorldMat` and :meth:`worldToVoxMat` attributes contain
-    transformation matrices for transforming between voxel and world
-    coordinates. The ``Nifti`` class follows the same process as ``nibabel``
-    in selecting the affine (see
+    The ``Nifti`` class is aware of three coordinate systems:
+
+      - The ``voxel`` coordinate system, used to access image data
+
+      - The ``world`` coordinate system, where voxel coordinates are
+        transformed into a millimetre coordinate system, defined by the
+        ``sform`` and/or ``qform`` elements of the NIFTI header.
+
+      - The ``fsl`` coordinate system, where voxel coordinated are scaled by
+        the ``pixdim`` values in the NIFTI header, and the X axis is inverted
+        if the voxel-to-world affine has a positive determinant.
+
+
+    The :meth:`getAffine` method is a simple means of acquiring an affine
+    which will transform between any of these coordinate systems.
+
+
+    See `here <http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the_format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to_the_transformation_parameters.3F>`_
+    for more details on the ``fsl`` coordinate system..
+
+
+    The ``Nifti`` class follows the same process as ``nibabel`` in determining
+    the ``voxel`` to ``world`` affine (see
     http://nipy.org/nibabel/nifti_images.html#the-nifti-affines):
 
 
@@ -227,20 +246,51 @@ class Nifti(notifier.Notifier, meta.Meta):
         header                   = header
         origShape, shape, pixdim = self.__determineShape(header)
 
-        voxToWorldMat = self.__determineTransform(header)
-        worldToVoxMat = transform.invert(voxToWorldMat)
+        self.header      = header
+        self.__shape     = shape
+        self.__origShape = origShape
+        self.__pixdim    = pixdim
+        self.__intent    = header.get('intent_code',
+                                      constants.NIFTI_INTENT_NONE)
+
+        voxToWorldMat    = self.__determineAffine(header)
+        affines, isneuro = self.__generateAffines(voxToWorldMat)
 
-        self.header          = header
-        self.__shape         = shape
-        self.__intent        = header.get('intent_code',
-                                          constants.NIFTI_INTENT_NONE)
-        self.__origShape     = origShape
-        self.__pixdim        = pixdim
-        self.__voxToWorldMat = voxToWorldMat
-        self.__worldToVoxMat = worldToVoxMat
+        self.__affines        = affines
+        self.__isNeurological = isneuro
 
 
-    def __determineTransform(self, header):
+    def __determineShape(self, header):
+        """This method is called by :meth:`__init__`. It figures out the actual
+        shape of the image data, and the zooms/pixdims for each data axis. Any
+        empty trailing dimensions are squeezed, but the returned shape is
+        guaranteed to be at least 3 dimensions. Returns:
+
+         - A sequence/tuple containing the image shape, as reported in the
+           header.
+         - A sequence/tuple containing the effective image shape.
+         - A sequence/tuple containing the zooms/pixdims.
+        """
+
+        # The canonicalShape function figures out
+        # the data shape that we should use.
+        origShape = list(header.get_data_shape())
+        shape     = canonicalShape(origShape)
+        pixdims   = list(header.get_zooms())
+
+        # if get_zooms() doesn't return at
+        # least len(shape) values, we'll
+        # fall back to the pixdim field in
+        # the header.
+        if len(pixdims) < len(shape):
+            pixdims = header['pixdim'][1:]
+
+        pixdims = pixdims[:len(shape)]
+
+        return origShape, shape, pixdims
+
+
+    def __determineAffine(self, header):
         """Called by :meth:`__init__`. Figures out the voxel-to-world
         coordinate transformation matrix that is associated with this
         ``Nifti`` instance.
@@ -250,26 +300,9 @@ class Nifti(notifier.Notifier, meta.Meta):
         # specially, as FNIRT clobbers the
         # sform section of the NIFTI header
         # to store other data.
-        #
-        # TODO Nibabel <= 2.1.0 has a bug in header.get
-        #      for fields with a value of 0. When this
-        #      bug gets fixed, we can replace the if-else
-        #      block below with this:
-        #
-        #          intent = header.get('intent_code', -1)
-        #          qform  = header.get('qform_code',  -1)
-        #          sform  = header.get('sform_code',  -1)
-        #
-        # TODO Change this in fslpy 2.0.0
-        #
-        if isinstance(header, nib.nifti1.Nifti1Header):
-            intent = header['intent_code']
-            qform  = header['qform_code']
-            sform  = header['sform_code']
-        else:
-            intent = -1
-            qform  = -1
-            sform  = -1
+        intent = header.get('intent_code', -1)
+        qform  = header.get('qform_code',  -1)
+        sform  = header.get('sform_code',  -1)
 
         if intent in (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
                       constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
@@ -300,34 +333,47 @@ class Nifti(notifier.Notifier, meta.Meta):
         return voxToWorldMat
 
 
-    def __determineShape(self, header):
-        """This method is called by :meth:`__init__`. It figures out the actual
-        shape of the image data, and the zooms/pixdims for each data axis. Any
-        empty trailing dimensions are squeezed, but the returned shape is
-        guaranteed to be at least 3 dimensions. Returns:
+    def __generateAffines(self, voxToWorldMat):
+        """Called by :meth:`__init__`, and the :meth:`voxToWorldMat` setter.
+        Generates and returns a dictionary containing affine transformations
+        between the ``voxel``, ``fsl``, and ``world`` coordinate
+        systems. These affines are accessible via the :meth:`getAffine`
+        method.
 
-         - A sequence/tuple containing the image shape, as reported in the
-           header.
-         - A sequence/tuple containing the effective image shape.
-         - A sequence/tuple containing the zooms/pixdims.
+        :arg voxToWorldMat: The voxel-to-world affine transformation
+        :returns:           A tuple containing:
+
+                             - a dictionary of affine transformations between
+                               each pair of coordinate systems
+                             - ``True`` if the image is to be considered
+                               "neurological", ``False`` otherwise - see the
+                               :meth:`isNeurological` method.
         """
 
-        # The canonicalShape function figures out
-        # the data shape that we should use.
-        origShape = list(header.get_data_shape())
-        shape     = canonicalShape(origShape)
-        pixdims   = list(header.get_zooms())
+        import numpy.linalg as npla
 
-        # if get_zooms() doesn't return at
-        # least len(shape) values, we'll
-        # fall back to the pixdim field in
-        # the header.
-        if len(pixdims) < len(shape):
-            pixdims = header['pixdim'][1:]
+        affines = {}
+        shape             = list(self.shape[ :3])
+        pixdim            = list(self.pixdim[:3])
+        voxToScaledVoxMat = np.diag(pixdim + [1.0])
+        isneuro           = npla.det(voxToWorldMat) > 0
 
-        pixdims = pixdims[:len(shape)]
+        if isneuro:
+            x                 = (shape[0] - 1) * pixdim[0]
+            flip              = transform.scaleOffsetXform([-1, 1, 1],
+                                                           [ x, 0, 0])
+            voxToScaledVoxMat = transform.concat(flip, voxToScaledVoxMat)
 
-        return origShape, shape, pixdims
+        affines['voxel', 'world'] = voxToWorldMat
+        affines['world', 'voxel'] = transform.invert(voxToWorldMat)
+        affines['voxel', 'fsl']   = voxToScaledVoxMat
+        affines['fsl',   'voxel'] = transform.invert(voxToScaledVoxMat)
+        affines['fsl',   'world'] = transform.concat(affines['voxel', 'world'],
+                                                     affines['fsl',   'voxel'])
+        affines['world', 'fsl']   = transform.concat(affines['fsl',   'voxel'],
+                                                     affines['voxel', 'world'])
+
+        return affines, isneuro
 
 
     def strval(self, key):
@@ -382,6 +428,16 @@ class Nifti(notifier.Notifier, meta.Meta):
         return tuple(self.__shape)
 
 
+    @property
+    def ndim(self):
+        """Returns the number of dimensions in this image. This number may not
+        match the number of dimensions specified in the NIFTI header, as
+        trailing dimensions of length 1 are ignored. But it is guaranteed to be
+        at least 3.
+        """
+        return len(self.__shape)
+
+
     @property
     def pixdim(self):
         """Returns a tuple containing the image pixdims (voxel sizes)."""
@@ -426,12 +482,39 @@ class Nifti(notifier.Notifier, meta.Meta):
         return units
 
 
+    def getAffine(self, from_, to):
+        """Return an affine transformation which can be used to transform
+        coordinates from ``from_`` to ``to``.
+
+        Valid values for the ``from_`` and ``to`` arguments are:
+
+         - ``'voxel'``: The voxel coordinate system
+         - ``'world'``: The world coordinate system, as defined by the image
+           sform/qform
+         - ``'fsl'``: The FSL coordinate system (scaled voxels, with a
+           left-right flip if the sform/qform has a positive determinant)
+
+        :arg from_: Source coordinate system
+        :arg to:    Destination coordinate system
+        :returns:   A ``numpy`` array of shape ``(4, 4)``
+        """
+        from_ = from_.lower()
+        to    = to   .lower()
+
+        if from_ not in ('voxel', 'fsl', 'world') or \
+           to    not in ('voxel', 'fsl', 'world'):
+            raise ValueError('Invalid source/reference spaces: '
+                             '{} -> {}'.format(from_, to))
+
+        return np.copy(self.__affines[from_, to])
+
+
     @property
     def worldToVoxMat(self):
         """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
         affine transformation from world coordinates to voxel coordinates.
         """
-        return np.array(self.__worldToVoxMat)
+        return self.getAffine('world', 'voxel')
 
 
     @property
@@ -439,7 +522,7 @@ class Nifti(notifier.Notifier, meta.Meta):
         """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
         affine transformation from voxel coordinates to world coordinates.
         """
-        return np.array(self.__voxToWorldMat)
+        return self.getAffine('voxel', 'world')
 
 
     @voxToWorldMat.setter
@@ -463,8 +546,10 @@ class Nifti(notifier.Notifier, meta.Meta):
 
         header.set_sform(xform, code=sformCode)
 
-        self.__voxToWorldMat = self.__determineTransform(header)
-        self.__worldToVoxMat = transform.invert(self.__voxToWorldMat)
+        affines, isneuro = self.__generateAffines(xform)
+
+        self.__affines        = affines
+        self.__isNeurological = isneuro
 
         log.debug('Affine changed:\npixdims: '
                   '%s\nsform: %s\nqform: %s',
@@ -475,6 +560,27 @@ class Nifti(notifier.Notifier, meta.Meta):
         self.notify(topic='transform')
 
 
+    @property
+    def voxToScaledVoxMat(self):
+        """Returns a transformation matrix which transforms from voxel
+        coordinates into scaled voxel coordinates, with a left-right flip
+        if the image appears to be stored in neurological order.
+
+        See http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the\
+        _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
+        _the_transformation_parameters.3F
+        """
+        return self.getAffine('voxel', 'fsl')
+
+
+    @property
+    def scaledVoxToVoxMat(self):
+        """Returns a transformation matrix which transforms from scaled voxels
+        into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform.
+        """
+        return self.getAffine('fsl', 'voxel')
+
+
     def mapIndices(self, sliceobj):
         """Adjusts the given slice object so that it may be used to index the
         underlying ``nibabel`` NIFTI image object.
@@ -490,16 +596,6 @@ class Nifti(notifier.Notifier, meta.Meta):
         return fileslice.canonical_slicers(sliceobj, self.__origShape)
 
 
-    @property
-    def ndim(self):
-        """Returns the number of dimensions in this image. This number may not
-        match the number of dimensions specified in the NIFTI header, as
-        trailing dimensions of length 1 are ignored. But it is guaranteed to be
-        at least 3.
-        """
-        return len(self.__shape)
-
-
     def getXFormCode(self, code=None):
         """This method returns the code contained in the NIFTI header,
         indicating the space to which the (transformed) image is oriented.
@@ -546,6 +642,7 @@ class Nifti(notifier.Notifier, meta.Meta):
 
         return int(code)
 
+
     # TODO Check what has worse performance - hashing
     #      a 4x4 array (via memoizeMD5), or the call
     #      to aff2axcodes. I'm guessing the latter,
@@ -563,7 +660,6 @@ class Nifti(notifier.Notifier, meta.Meta):
         return nib.orientations.aff2axcodes(xform, inaxes)
 
 
-    @memoize.Instanceify(memoize.memoize)
     def isNeurological(self):
         """Returns ``True`` if it looks like this ``Nifti`` object has a
         neurological voxel orientation, ``False`` otherwise. This test is
@@ -582,51 +678,7 @@ class Nifti(notifier.Notifier, meta.Meta):
         _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
         _the_transformation_parameters.3F
         """
-        import numpy.linalg as npla
-        return npla.det(self.__voxToWorldMat) > 0
-
-
-    @property
-    def voxToScaledVoxMat(self):
-        """Returns a transformation matrix which transforms from voxel
-        coordinates into scaled voxel coordinates, with a left-right flip
-        if the image appears to be stored in neurological order.
-
-        See http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FLIRT/FAQ#What_is_the\
-        _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\
-        _the_transformation_parameters.3F
-        """
-        return self.__voxToScaledVoxMat()
-
-
-    @memoize.Instanceify(memoize.memoize)
-    def __voxToScaledVoxMat(self):
-        """See :meth:`voxToScaledVoxMat`. """
-
-        shape          = list(self.shape[ :3])
-        pixdim         = list(self.pixdim[:3])
-        voxToPixdimMat = np.diag(pixdim + [1.0])
-
-        if self.isNeurological():
-            x              = (shape[0] - 1) * pixdim[0]
-            flip           = transform.scaleOffsetXform([-1, 1, 1], [x, 0, 0])
-            voxToPixdimMat = transform.concat(flip, voxToPixdimMat)
-
-        return voxToPixdimMat
-
-
-    @property
-    def scaledVoxToVoxMat(self):
-        """Returns a transformation matrix which transforms from scaled voxels
-        into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform.
-        """
-        return self.__scaledVoxToVoxMat()
-
-
-    @memoize.Instanceify(memoize.memoize)
-    def __scaledVoxToVoxMat(self):
-        """See :meth:`scaledVoxToVoxMat`. """
-        return transform.invert(self.voxToScaledVoxMat)
+        return self.__isNeurological
 
 
     def sameSpace(self, other):