From 71843960e2b6ced19c74e0a8356dd566429d151e Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauld.mccarthy@gmail.com>
Date: Thu, 12 Jan 2017 15:49:24 +0000
Subject: [PATCH] Nifti instances now allow voxToWorldMat to be edited, and
notify listeners when it is. New functions in transform module.
---
fsl/data/image.py | 171 +++++++++++++++++++++++++++++++----------
fsl/utils/transform.py | 104 +++++++++++++++++++++++++
2 files changed, 233 insertions(+), 42 deletions(-)
diff --git a/fsl/data/image.py b/fsl/data/image.py
index eae93b74e..86b16d908 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -51,7 +51,7 @@ import fsl.data.imagewrapper as imagewrapper
log = logging.getLogger(__name__)
-class Nifti(object):
+class Nifti(notifier.Notifier):
"""The ``Nifti`` class is intended to be used as a base class for
things which either are, or are associated with, a NIFTI image.
The ``Nifti`` class is intended to represent information stored in
@@ -104,6 +104,15 @@ class Nifti(object):
:attr:`.constants.NIFTI_XFORM_ANALYZE`.
+ The ``Nifti`` class implements the :class:`.Notifier` interface -
+ listeners may register to be notified on the following topics:
+
+ =============== ========================================================
+ ``'transform'`` The affine transformation matrix has changed. This topic
+ will occur when the ``voxToWorldMat`` is changed.
+ =============== ========================================================
+
+
.. warning:: The ``header`` field may either be a ``nifti1``, ``nifti2``,
or ``analyze`` header object. Make sure to take this into
account if you are writing code that should work with all
@@ -144,37 +153,16 @@ class Nifti(object):
voxToWorldMat = self.__determineTransform(header)
worldToVoxMat = transform.invert(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
-
-
- @property
- def niftiVersion(self):
- """Returns the NIFTI file version:
-
- - ``0`` for ANALYZE
- - ``1`` for NIFTI1
- - ``2`` for NIFTI2
- """
-
- import nibabel as nib
+ 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
- # nib.Nifti2 is a subclass of Nifti1,
- # and Nifti1 a subclass of Analyze,
- # so we have to check in order
- if isinstance(self.header, nib.nifti2.Nifti2Header): return 2
- elif isinstance(self.header, nib.nifti1.Nifti1Header): return 1
- elif isinstance(self.header, nib.analyze.AnalyzeHeader): return 0
-
- else: raise RuntimeError('Unrecognised header: {}'.format(self.header))
-
-
+
def __determineTransform(self, header):
"""Called by :meth:`__init__`. Figures out the voxel-to-world
coordinate transformation matrix that is associated with this
@@ -248,6 +236,87 @@ class Nifti(object):
return origShape, shape, pixdims
+ @property
+ def niftiVersion(self):
+ """Returns the NIFTI file version:
+
+ - ``0`` for ANALYZE
+ - ``1`` for NIFTI1
+ - ``2`` for NIFTI2
+ """
+
+ import nibabel as nib
+
+ # nib.Nifti2 is a subclass of Nifti1,
+ # and Nifti1 a subclass of Analyze,
+ # so we have to check in order
+ if isinstance(self.header, nib.nifti2.Nifti2Header): return 2
+ elif isinstance(self.header, nib.nifti1.Nifti1Header): return 1
+ elif isinstance(self.header, nib.analyze.AnalyzeHeader): return 0
+
+ else: raise RuntimeError('Unrecognised header: {}'.format(self.header))
+
+
+ @property
+ def shape(self):
+ """Returns a tuple containing the image data shape. """
+ return tuple(self.__shape)
+
+
+ @property
+ def pixdim(self):
+ """Returns a tuple containing the image pixdims (voxel sizes)."""
+ return tuple(self.__pixdim)
+
+
+ @property
+ def intent(self):
+ """Returns the NIFTI intent code of this image.
+ """
+ return self.__intent
+
+
+ @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)
+
+
+ @property
+ def voxToWorldMat(self):
+ """Returns a ``numpy`` array of shape ``(4, 4)`` containing an
+ affine transformation from voxel coordinates to world coordinates.
+ """
+ return np.array(self.__voxToWorldMat)
+
+
+ @voxToWorldMat.setter
+ def voxToWorldMat(self, xform):
+ """Update the ``voxToWorldMat``. The ``worldToVoxMat`` and ``pixdim``
+ values are also updated. This will result in notification on the
+ ``'transform'`` topic.
+ """
+
+ header = self.header
+
+ header.set_qform(xform)
+ header.set_sform(xform)
+
+ self.__voxToWorldMat = self.__determineTransform(header)
+ self.__worldToVoxMat = transform.invert(self.__voxToWorldMat)
+ self.__pixdim = header.get_zooms()
+
+ log.debug('{} affine changed:\npixdims: '
+ '{}\nsform: {}\nqform: {}'.format(
+ header.get_zooms(),
+ header.get_sform(),
+ header.get_qform()))
+
+ self.notify(topic='transform')
+
+
def mapIndices(self, sliceobj):
"""Adjusts the given slice object so that it may be used to index the
underlying ``nibabel`` NIFTI image object.
@@ -267,7 +336,7 @@ class Nifti(object):
# TODO: Remove this method, and use the shape attribute directly
def is4DImage(self):
"""Returns ``True`` if this image is 4D, ``False`` otherwise. """
- return len(self.shape) > 3 and self.shape[3] > 1
+ return len(self.__shape) > 3 and self.__shape[3] > 1
def getXFormCode(self, code=None):
@@ -350,7 +419,7 @@ class Nifti(object):
_the_transformation_parameters.3F
"""
import numpy.linalg as npla
- return npla.det(self.voxToWorldMat) > 0
+ return npla.det(self.__voxToWorldMat) > 0
def sameSpace(self, other):
@@ -358,9 +427,12 @@ class Nifti(object):
:class:`Nifti` instance) has the same dimensions and is in the
same space as this image.
"""
- return np.all(np.isclose(self.shape[:3], other.shape[:3])) and \
- np.all(np.isclose(self.pixdim, other.pixdim)) and \
- np.all(np.isclose(self.voxToWorldMat, other.voxToWorldMat))
+ return np.all(np.isclose(self .__shape[:3],
+ other.__shape[:3])) and \
+ np.all(np.isclose(self .__pixdim,
+ other.__pixdim)) and \
+ np.all(np.isclose(self .__voxToWorldMat,
+ other.__voxToWorldMat))
def getOrientation(self, axis, xform):
@@ -405,7 +477,7 @@ class Nifti(object):
return code
-class Image(Nifti, notifier.Notifier):
+class Image(Nifti):
"""Class which represents a 3D/4D NIFTI image. Internally, the image is
loaded/stored using a :mod:`nibabel.nifti1.Nifti1Image` or
:mod:`nibabel.nifti2.Nifti2Image`, and data access managed by a
@@ -439,10 +511,10 @@ class Image(Nifti, notifier.Notifier):
============== ===========================================================
- The ``Image`` class implements the :class:`.Notifier` interface -
- listeners may register to be notified of changes to the above properties,
- by registering on the following _topic_ names (see the :class:`.Notifier`
- class documentation):
+ The ``Image`` class adds some :class:`.Notifier` topics to those which are
+ already provided by the :class:`Nifti` class - listeners may register to
+ be notified of changes to the above properties, by registering on the
+ following _topic_ names (see the :class:`.Notifier` class documentation):
=============== ======================================================
@@ -453,7 +525,8 @@ class Image(Nifti, notifier.Notifier):
value (see :meth:`.Notifier.notify`).
``'saveState'`` This topic is notified whenever the saved state of the
- image changes (i.e. it is edited, or saved to disk).
+ image changes (i.e. data or ``voxToWorldMat`` is
+ edited, or the image saved to disk).
``'dataRange'`` This topic is notified whenever the image data range
is changed/adjusted.
@@ -596,6 +669,11 @@ class Image(Nifti, notifier.Notifier):
loadData=loadData,
threaded=threaded)
+ # Listen to ourself for changes
+ # to the voxToWorldMat, so we
+ # can update the saveState.
+ self.register(self.name, self.__transformChanged, topic='transform')
+
if calcRange:
self.calcRange()
@@ -688,6 +766,15 @@ class Image(Nifti, notifier.Notifier):
return self.__nibImage.dataobj[tuple(coords)].dtype
+ def __transformChanged(self, *args, **kwargs):
+ """Called when the ``voxToWorldMat`` of this :class:`Nifti` instance
+ changes. Updates the :attr:`saveState` accordinbgly.
+ """
+ if self.__saveState:
+ self.__saveState = False
+ self.notify(topic='saveState')
+
+
def __dataRangeChanged(self, *args, **kwargs):
"""Called when the :class:`.ImageWrapper` data range changes.
Notifies any listeners of this ``Image`` (registered through the
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index bb7a96af5..ed90bf35e 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -14,6 +14,10 @@ spaces. The following functions are provided:
scaleOffsetXform
invert
concat
+ compse
+ decompose
+ rotMatToAxisAngles
+ axisAnglesToRotMat
axisBounds
"""
@@ -76,6 +80,106 @@ def scaleOffsetXform(scales, offsets):
return xform
+def compose(scales, offsets, rotations, origin=None):
+ """Compose a transformation matrix out of the given scales, offsets
+ and axis rotations.
+
+ :arg scales: Sequence of three scale values.
+
+ :arg offsets: Sequence of three offset values.
+
+ :arg rotations: Sequence of three rotation values, in radians.
+
+ :arg origin: Origin of rotation - must be scaled by the ``scales``.
+ If not provided, the rotation origin is ``(0, 0, 0)``.
+ """
+
+ preRotate = np.eye(4)
+ postRotate = np.eye(4)
+ if origin is not None:
+ preRotate[ 0, 3] = -origin[0]
+ preRotate[ 1, 3] = -origin[1]
+ preRotate[ 2, 3] = -origin[2]
+ postRotate[0, 3] = origin[0]
+ postRotate[1, 3] = origin[1]
+ postRotate[2, 3] = origin[2]
+
+ scale = np.eye(4, dtype=np.float64)
+ offset = np.eye(4, dtype=np.float64)
+ rotate = np.eye(4, dtype=np.float64)
+
+ scale[ 0, 0] = scales[ 0]
+ scale[ 1, 1] = scales[ 1]
+ scale[ 2, 2] = scales[ 2]
+ offset[ 0, 3] = offsets[0]
+ offset[ 1, 3] = offsets[1]
+ offset[ 2, 3] = offsets[2]
+ rotate[:3, :3] = axisAnglesToRotMat(*rotations)
+
+ return concat(offset, postRotate, rotate, preRotate, scale)
+
+
+def decompose(xform):
+ """Decomposes the given transformation matrix into separate offsets,
+ scales, and rotations.
+
+ .. note:: Shears are not yet supported, and may never be supported.
+ """
+
+ offsets = xform[:3, 3]
+ scales = [np.sqrt(np.sum(xform[:3, 0] ** 2)),
+ np.sqrt(np.sum(xform[:3, 1] ** 2)),
+ np.sqrt(np.sum(xform[:3, 2] ** 2))]
+
+ rotmat = np.copy(xform[:3, :3])
+ rotmat[:, 0] /= scales[0]
+ rotmat[:, 1] /= scales[1]
+ rotmat[:, 2] /= scales[2]
+
+ rots = rotMatToAxisAngles(rotmat)
+
+ return scales, offsets, rots
+
+
+def rotMatToAxisAngles(rotmat):
+ """Given a ``(3, 3)`` rotation matrix, decomposes the rotations into
+ an angle in radians about each axis.
+ """
+ xrot = np.arctan2(rotmat[2, 1], rotmat[2, 2])
+ yrot = np.sqrt( rotmat[2, 1] ** 2 + rotmat[2, 2] ** 2)
+ yrot = np.arctan2(rotmat[2, 0], yrot)
+ zrot = np.arctan2(rotmat[1, 0], rotmat[0, 0])
+
+ return [xrot, yrot, zrot]
+
+
+def axisAnglesToRotMat(xrot, yrot, zrot):
+ """Constructs a ``(3, 3)`` rotation matrix from the given angles, which
+ must be specified in radians.
+ """
+
+ xmat = np.eye(3)
+ ymat = np.eye(3)
+ zmat = np.eye(3)
+
+ xmat[1, 1] = np.cos(xrot)
+ xmat[1, 2] = -np.sin(xrot)
+ xmat[2, 1] = np.sin(xrot)
+ xmat[2, 2] = np.cos(xrot)
+
+ ymat[0, 0] = np.cos(yrot)
+ ymat[0, 2] = np.sin(yrot)
+ ymat[2, 0] = -np.sin(yrot)
+ ymat[2, 2] = np.cos(yrot)
+
+ zmat[0, 0] = np.cos(zrot)
+ zmat[0, 1] = -np.sin(zrot)
+ zmat[1, 0] = np.sin(zrot)
+ zmat[1, 1] = np.cos(zrot)
+
+ return concat(zmat, ymat, xmat)
+
+
def axisBounds(shape,
xform,
axes=None,
--
GitLab