From 0b1b35bb7f31df4bf4043112f2f19d7bdb063f6e Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Tue, 9 Apr 2019 12:13:32 +0100
Subject: [PATCH] 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 | 288 +++++++++++++++++++++++++++-------------------
1 file changed, 170 insertions(+), 118 deletions(-)
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 1bf901c28..6832e560f 100644
--- a/fsl/data/image.py
+++ b/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):
--
GitLab