diff --git a/fsl/transform/x5.py b/fsl/transform/x5.py
index b37dee2d75f6eeeec3049b2c9447cf7140d81e03..b02894f0fa245ca313c956d666eecaceaa7595e3 100644
--- a/fsl/transform/x5.py
+++ b/fsl/transform/x5.py
@@ -5,11 +5,19 @@
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains functions for reading/writing linear/non-linear FSL
-transformations from/to BIDS X5 files.
+transformations from/to BIDS X5 files. The following functions are available:
+
+.. autosummary::
+ :nosignatures:
+
+ readLinearX5
+ writeLinearX5
+ readNonLinearX5
+ writeNonLinearX5
An X5 file is a HDF5 container file which stores a linear or non-linear
-transformation between a **source** NIfTI image and a **refernece** image. In
+transformation between a **source** NIfTI image and a **reference** image. In
an X5 file, the source image is referred to as the ``From`` space, and the
reference image the ``To`` space.
@@ -18,13 +26,16 @@ Custom HDF5 groups
------------------
-HDF5 files are composed of *groups*, *attributes*, and *datasets*. Groups are
-simply containers for attributes, datasets, and other groups. Attributes are
-strongly-typed scalar values. Datasets are strongly-typed, structured
-N-dimensional arrays.
+HDF5 files are composed primarily of *groups*, *attributes*, and
+*datasets*:
+
+ - *Groups* are simply containers for attributes, datasets, and other groups.
+ - *Attributes* are strongly-typed scalar values.
+ - *Datasets* are strongly-typed, structured N-dimensional arrays.
-To simplify the file format efinitions below, we shall first define a few
+
+To simplify the file format definitions below, we shall first define a few
custom HDF5 groups. In the file format definitions, a HDF5 group which is
listed as being of one of these custom types shall contain the
attributes/datasets that are listed here.
@@ -40,37 +51,43 @@ another. Groups of type "affine" have the following fields:
+---------------+-----------+-------------------------------------------------+
-| Name | Type | Value/Description |
+| **Name** | **Type** | **Value/Description** |
+---------------+-----------+-------------------------------------------------+
| ``Type`` | attribute | ``'linear'`` |
-| ``Transform`` | dataset | The transformation itself - a 4x4 ``float64`` |
-| | | affine transformation |
++---------------+-----------+-------------------------------------------------+
+| ``Transform`` | dataset | The affine transformation - a ``float64`` array |
+| | | of shape ``(4, 4)`` |
++---------------+-----------+-------------------------------------------------+
| ``Inverse`` | dataset | Optional pre-calculated inverse |
+---------------+-----------+-------------------------------------------------+
-*mapping*
-^^^^^^^^^
+*space*
+^^^^^^^
-A HDF5 group which is listed as being of type "mapping" contains all of the
+A HDF5 group which is listed as being of type "space" contains all of the
information required to define the space of a NIfTI image, including its
-dimensions, and voxel-to-world affine transformation. The *world* coordinate
-system of an image is defined by its ``sform`` or ``qform`` (hereafter
-referred to as the ``sform``), which is contained in the NIfTI header.
+shape, dimensions, and voxel-to-world affine transformation. The *world*
+coordinate system of an image is defined by its ``sform`` or ``qform``
+(hereafter referred to as the ``sform``), which is contained in the NIfTI
+header.
-Groups of type "mapping" have the following fields:
+Groups of type "space" have the following fields:
+-------------+-----------+---------------------------------------------------+
-| Name | Type | Value/Description |
+| **Name** | **Type** | **Value/Description** |
+-------------+-----------+---------------------------------------------------+
| ``Type`` | attribute | ``'image'`` |
++-------------+-----------+---------------------------------------------------+
| ``Size`` | attribute | ``uint64`` ``(X, Y, Z)`` voxel dimensions |
++-------------+-----------+---------------------------------------------------+
| ``Scales`` | attribute | ``float64`` ``(X, Y, Z)`` voxel pixdims |
++-------------+-----------+---------------------------------------------------+
| ``Mapping`` | affine | The image voxel-to-world transformation (its |
-| | | "sform") |
+| | | ``sform``) |
+-------------+-----------+---------------------------------------------------+
@@ -80,7 +97,7 @@ Linear X5 files
Linear X5 transformation files contain an affine transformation matrix of
shape ``(4, 4)``, which can be used to transform **source image world
-coordinates into reference image world coordinates**.
+coordinates** into **reference image world coordinates**.
Linear X5 transformation files are assumed to adhere to the HDF5 structure
@@ -88,23 +105,27 @@ defined in the table below. All fields are required.
+-----------------------------+-----------+-----------------------------------+
-| Name | Type | Value/Description |
+| **Name** | **Type** | **Value/Description** |
+-----------------------------+-----------+-----------------------------------+
| *Metadata* |
+-----------------------------+-----------+-----------------------------------+
| ``/Format`` | attribute | ``'X5'`` |
++-----------------------------+-----------+-----------------------------------+
| ``/Version`` | attribute | ``'0.0.1'`` |
++-----------------------------+-----------+-----------------------------------+
| ``/Metadata`` | attribute | JSON string containing |
| | | unstructured metadata. |
+-----------------------------+-----------+-----------------------------------+
| *Transformation* |
+-----------------------------+-----------+-----------------------------------+
-| ``/From/`` | mapping | Source image mapping |
-| ``/To/`` | mapping | Reference image mapping |
| ``/`` | affine | Affine transformation from source |
| | | image world coordinates to |
| | | reference image world coordinates |
+-----------------------------+-----------+-----------------------------------+
+| ``/From/`` | space | Source image definition |
++-----------------------------+-----------+-----------------------------------+
+| ``/To/`` | space | Reference image definition |
++-----------------------------+-----------+-----------------------------------+
Non-linear X5 transformation files
@@ -116,71 +137,195 @@ a source image coordinate system and a reference image coordinate system. The
transformation is represented as either:
- A *displacement field*, which is defined in the same space as the reference
- image, and which contains displacements.
+ image, and which contains displacements from reference image coordinates to
+ source image coordinates.
- A quadratic or cubic B-spline *coefficient field*, which contains
- coefficients defined in a coarse grid which overlays the same space as the
+ coefficients defined in a coarse grid overlaid onto the same space as the
reference image, and from which a displacement field can be calculated.
-Displacement fields may contain either:
- - *relative displacements*, which contains displacements *from* the
- reference image coordinates *to* the source image coordinates (i.e. add the
- displacements to the reference image coordinates to get the corresponding
- source image coordinates)
+Displacement fields
+^^^^^^^^^^^^^^^^^^^
+
+
+A displacement field is a ``float64`` array of shape ``(X, Y, Z, 3)``, defined
+in the same space as the reference image (i.e. the reference image is assumed
+to have shape ``(X, Y, Z)``. A displacement field may contain either:
+
+ - *relative* displacements, where each voxel in the displacement field
+ contains an offset which can be added to the reference image coordinates
+ for that voxel, in order to calculate the corresponding source image
+ coordinates.
+
+ - *absolute* displacements, where each voxel in the displacement field simply
+ contains the source image coordinates which correspond to those reference
+ image coordinates.
+
+
+Coefficient fields
+^^^^^^^^^^^^^^^^^^
+
+
+A coefficient field is a ``float64`` array of shape ``(X, Y, Z, 3)`` which
+contains the coefficients of a set of quadratic or cubic B-spline functions
+defined on a regular 3D grid overlaid on the reference image voxel coordinate
+system. Each coefficient in this grid may be referred to as a *control point*
+or a *knot*.
+
+
+Evaluating the spline functions at a particular location in the grid will
+result in a relative displacement which can be applied to that location's
+reference image coordinates, in order to determine the corresponding source
+image coordinates.
+
+
+The shape of this coefficient grid is not necessarily the same as the shape of
+the reference image grid. For this reason, some additional parameters are
+stored in coefficient field files, in a sub-group called ``/Parameters``:
+
+ - The distance between control points, defined in terms of reference image
+ voxels
+ - An affine transformation which can be used to transform reference image
+ voxel coordinates into coefficient field voxel coordinates.
+
+
+Non-linear coordinate systems
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+The coordinate systems used in a displacement field, or in a displacement
+field that has been generated from a coefficient field, defines relative
+displacements from a reference image coordinate system to a source image
+coordinate system. The coordinate systems used are not defined in this
+specification - they may be voxels, world coordinates, FSL coordinates, or
+some other coordinate system.
+
+
+Howewer, if the transformation does not transform between source and reference
+image **world** coordinates, the ``/Pre`` and ``/Post`` affine transformations
+must be provided.
+
+The ``/Pre/`` affine transformation will be used to transform reference image
+world coordinates into the reference image coordinate system required for use
+with the displacement values (or required for the spline coefficient
+evaluation). Similarly, the ``/Post/`` affine transformation will be used to
+transform the resulting source image coordinates into the source image world
+coordinate system.
+
+
+Initial affine alignment
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+Non-linear transformations are often accompanied by an initial affine
+transformation, which provides a coarse global initial alignment that is used
+as the starting point for the non-linear registration process. The non-linear
+transformation will then typically encode displacements between the
+registered, or aligned source image, and the reference image.
- - *absolute displacement* which simply contain the source image coordinates
- at each voxel in the reference image space.
+Now we have three spaces, and three sets of coordinate systems, to consider:
-Displacement field transformations define displacements between the refernece
-image *world* coordinate system, and the source image *world* coordinate
-system.
+ 1. Source image space - the source image, before initial linear registration
+ to the reference image
+ 2. Aligned-source image space - the source image, after it has been linearly
+ transformed to the reference image space
-The displacement field generated from a coefficient field defines relative
-displacements from a reference image space to a source image space. These
-spaces are undefined; however, the ``/Pre`` affine transformation may be used
-to transform reference image world coordinates into the reference image space,
-and the ``/Post`` affine transformation may be used to transform the resulting
-source image coordinates into the source image world coordinate system.
+ 3. Reference image space
+
+The initial affine registration encodes a transformation between spaces 1 and
+2, and the non-linear registration encodes a transformation between spaces 2
+and 3. Note that the fields-of-view for spaces 2 and 3 are typically
+equivalent.
+
+
+The initial affine transformation may be included in an X5 file, in the
+``/InitialAlignment/`` group. If provided, this initial transformation is
+assumed to provide a transformation:
+
+ - *From* the **source image** world coordinate system (or the coordinate
+ system used as input to the ``/Post/`` affine, if provided).
+
+ - *To* the **aligned source** image coordinate system used within the
+ non-linear transformation.
+
+
+File format specification
+^^^^^^^^^^^^^^^^^^^^^^^^^
Non-linear X5 transformation files are assumed to adhere to the following
-HDF5 structure::
+HDF5 structure:
-+-----------------------------+-----------+-----------------------------------+
-| Name | Type | Value/Description |
-+-----------------------------+-----------+-----------------------------------+
-| *Metadata* - same as for linear transformation files |
-+-----------------------------------------------------------------------------+
++---------------------+-----------+-------------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++---------------------+-----------+-------------------------------------------+
+| *Metadata* |
++---------------------+-----------+-------------------------------------------+
+| ``/Format`` | attribute | ``'X5'`` |
++---------------------+-----------+-------------------------------------------+
+| ``/Version`` | attribute | ``'0.0.1'`` |
++---------------------+-----------+-------------------------------------------+
+| ``/Metadata`` | attribute | JSON string containing unstructured |
+| | | metadata. |
++---------------------+-----------+-------------------------------------------+
| *Transformation* |
-+-----------------------------+-----------+-----------------------------------+
-| ``/Type`` | attribute | ``'nonlinear'`` |
-| ``/SubType`` | attribute | ``'displacement'`` or |
-| | | ``'coefficient'`` |
-| ``/Representation`` | attribute | If ``/SubType`` is |
-| | | ``'displacement'``, its |
-| | | ``/Representation`` may be either |
-| | | ``'absolute'`` or ``'relative'``. |
-| | | If ``/SubType`` is |
-| | | ``'coefficient'``, its |
-| | | ``/Representation`` may be either |
-| | | ``'quadratic bspline'`` or |
-| | | ``'cubic bspline'``. |
-| ``/Transform`` | dataset | The transformation itself - see |
-| | | above. |
-| ``/Inverse`` | dataset | Optional pre-calculated inverse |
-| ``/From/`` | mapping | Source image mapping |
-| ``/To/`` | mapping | Reference image mapping |
-|-----------------------------+-----------+-----------------------------------+
++------------------------+-----------+----------------------------------------+
+| **Name** | **Type** | **Value/Description** |
++------------------------+-----------+----------------------------------------+
+| ``/Type`` | attribute | ``'nonlinear'`` |
++------------------------+-----------+----------------------------------------+
+| ``/SubType`` | attribute | ``'displacement'`` or |
+| | | ``'coefficient'`` |
++------------------------+-----------+----------------------------------------+
+| ``/Representation`` | attribute | If ``/SubType`` is ``'displacement'``, |
+| | | ``/Representation`` may be either |
+| | | ``'absolute'`` or ``'relative'``. |
+| | | If ``/SubType`` is ``'coefficient'``, |
+| | | ``/Representation`` may be either |
+| | | ``'quadratic bspline'`` or |
+| | | ``'cubic bspline'``. |
++------------------------+-----------+----------------------------------------+
+| ``/Transform`` | dataset | The displacement/coefficient field - |
+| | | see above description. |
++------------------------+-----------+----------------------------------------+
+| ``/Inverse`` | dataset | Optional pre-calculated inverse |
++------------------------+-----------+----------------------------------------+
+| ``/From/`` | space | Source image definition |
++------------------------+-----------+----------------------------------------+
+| ``/To/`` | space | Reference image definition |
++------------------------+-----------+----------------------------------------+
+| ``/Pre/`` | affine | Optional affine transformation from |
+| | | reference image world coordinates to |
+| | | the reference image coordinate system |
+| | | required by the displacement/ |
+| | | coefficient field. |
++------------------------+-----------+----------------------------------------+
+| ``/Post/`` | affine | Optional affine transformation from |
+| | | the source image coordinate system |
+| | | produced by the displacement/ |
+| | | coefficient field to the source image |
+| | | world coordinate system. |
++------------------------+-----------+----------------------------------------+
+| ``/InitialAlignment/`` | affine | Optional initial affine registration |
+| | | from the source image to the reference |
+| | | image. |
++------------------------+-----------+----------------------------------------+
| *Coefficient field parameters* (required for ``'coefficient'`` files) |
-|-----------------------------+-----------+-----------------------------------+
-| ``/Parameters/Spacing`` | attribute |
-| ``/Parameters/ReferenceToField`` | affine | | attribute |
-+-----------------------------+-----------+-----------------------------------+
-
++-----------------------------------+-----------+-----------------------------+
+| **Name** | **Type** | **Value/Description** |
++-----------------------------------+-----------+-----------------------------+
+| ``/Parameters/Spacing`` | attribute | ``uint64`` ``(X, Y, Z)`` |
+| | | knot spacing (defined in |
+| | | reference image voxels) |
++-----------------------------------+-----------+-----------------------------+
+| ``/Parameters/ReferenceToField/`` | affine | Reference image voxel to |
+| | | coefficient field voxel |
+| | | affine transformation. |
++-----------------------------------+-----------+-----------------------------+
"""
@@ -192,12 +337,18 @@ import h5py
import fsl.version as version
from . import affine
+from . import nonlinear
X5_FORMAT = 'X5'
X5_VERSION = '0.0.1'
+class InvalidFileError(Exception):
+ """Error raised if an invalid/incompatible file is detected. """
+ pass
+
+
def _writeMetadata(group):
"""Writes a metadata block into the given group.
@@ -208,8 +359,26 @@ def _writeMetadata(group):
group.attrs['Metadata'] = json.dumps({'fslpy' : version.__version__})
-def _readLinearTransform(group):
- """Reads a linear transformation from the given group,
+def _validateMetadata(group):
+ """Reads a metadata block from the given group, and raises an
+ :exc:`InvalidFileError` if it does not look valid.
+ """
+
+ try:
+ format = group.attrs['Format']
+ version = group.attrs['Version']
+
+ except Exception:
+ raise InvalidFileError('File does not appear to be an X5 file')
+
+ if (format != X5_FORMAT) or (version != X5_VERSION):
+ raise InvalidFileError('Incompatible format/version (required: {}/{}, '
+ 'present: {}/{})'.format(X5_FORMAT, X5_VERSION,
+ format, version))
+
+
+def _readAffine(group):
+ """Reads an affine from the given group,
:arg group: A ``h5py.Group`` object.
:returns: ``numpy`` array containing a ``(4, 4)`` affine transformation.
@@ -226,8 +395,8 @@ def _readLinearTransform(group):
return xform
-def _writeLinearTransform(group, xform):
- """Writes the given linear transformation and its inverse to the given
+def _writeAffine(group, xform):
+ """Writes the given affine transformation and its inverse to the given
group.
:arg group: A ``h5py.Group`` object.
@@ -242,9 +411,8 @@ def _writeLinearTransform(group, xform):
group.create_dataset('Inverse', data=inv)
-def _readLinearMapping(group):
- """Reads a linear mapping, defining a source or reference space, from the
- given group.
+def _readSpace(group):
+ """Reads a "space" group, defining a source or reference space.
:arg group: A ``h5py.Group`` object.
:returns: :class:`.Nifti` object. defining the mapping
@@ -257,7 +425,7 @@ def _readLinearMapping(group):
shape = group.attrs['Size']
pixdim = group.attrs['Scales']
- xform = _readLinearTransform(group['Mapping'])
+ xform = _readAffine(group['Mapping'])
hdr = nib.Nifti2Header()
hdr.set_data_shape(shape)
@@ -266,8 +434,8 @@ def _readLinearMapping(group):
return fslimage.Nifti(hdr)
-def _writeLinearMapping(group, img):
- """Writes a linear mapping specified by ``img`` to the given group.
+def _writeSpace(group, img):
+ """Writes a space specified by ``img`` to the given group.
:arg group: A ``h5py.Group`` object.
:arg img: :class:`.Nifti` object. defining the mapping
@@ -277,46 +445,100 @@ def _writeLinearMapping(group, img):
group.attrs['Scales'] = np.asarray(img.pixdim[:3], np.float32)
mapping = group.create_group('Mapping')
- _writeLinearTransform(mapping, img.getAffine('voxel', 'world'))
+ _writeAffine(mapping, img.getAffine('voxel', 'world'))
+
+def _readDisplacementField(group):
+ """
+ """
-def _readNonLinearTransform(group):
- if group.attrs['Type'] != 'nonlinear':
- raise ValueError('Not a nonlinear transform')
- return np.array(group['Transform'])
+ field = np.array(group['Transform'])
+ dispType = group.attrs['Representation']
+ src = _readSpace( group['From'])
+ ref = _readSpace( group['To'])
+ pre = _readAffine(group['Pre'])
+ post = _readAffine(group['Post'])
+ srcToRef = _readAffine(group['InitialAlignment'])
+ refSpace = affine.identify(ref, pre, from_='world')[1]
+ srcSpace = affine.identify(src, post, to='world')[ 0]
-def _writeNonLinearTransform(group, field):
+ return nonlinear.DisplacementField(field,
+ src,
+ ref,
+ srcSpace=srcSpace,
+ refSpace=refSpace,
+ srcToRefMat=srcToRef,
+ dispType=dispType)
+
+
+def _readCoefficientField(group):
"""
"""
- group.attrs['Type'] = 'nonlinear'
- group.create_dataset('Transform', data=field, dtype=np.float32)
+ field = np.array(group['Transform'])
+ fieldType = group.attrs['Representation']
+ src = _readSpace( group['From'])
+ ref = _readSpace( group['To'])
+ pre = _readAffine(group['Pre'])
+ post = _readAffine(group['Post'])
+ srcToRef = _readAffine(group['InitialAlignment'])
+
+ spacing = group['Parameters'].attrs['Spacing']
+ refToField = _readAffine(group['Parameters']['ReferenceToField'])
+ fieldToRef = affine.inverse(refToField)
+
+ refSpace = affine.identify(ref, pre, from_='world')[1]
+ srcSpace = affine.identify(src, post, to='world')[ 0]
+
+ return nonlinear.CoefficientField(field,
+ src,
+ ref,
+ srcSpace=srcSpace,
+ refSpace=refSpace,
+ fieldType=fieldType,
+ knotSpacing=spacing,
+ srcToRefMat=srcToRef,
+ fieldToRefMat=fieldToRef)
+
def readLinearX5(fname):
- """
+ """Read a linear X5 transformation file from ``fname``.
+
+ :return: A tuple containing:
+
+ - A ``(4, 4)`` ``numpy`` array containing the affine
+ transformation
+ - A :class:`.Nifti` instance representing the source space
+ - A :class:`.Nifti` instance representing the reference space
"""
with h5py.File(fname, 'r') as f:
- xform = _readLinearTransform(f['/'])
- src = _readLinearMapping( f['/From'])
- ref = _readLinearMapping( f['/To'])
+ _validateMetadata( f['/'])
+ xform = _readAffine( f['/'])
+ src = _readSpace( f['/From'])
+ ref = _readSpace( f['/To'])
return xform, src, ref
def writeLinearX5(fname, xform, src, ref):
- """
+ """Write a linear transformation to the specified file.
+
+ :arg fname: File name to write to
+ :arg xform: ``(4, 4)`` ``numpy`` array containing the affine transformation
+ :arg src: :class:`.Nifti` image representing the source space
+ :arg ref: :class:`.Nifti` image representing the reference
"""
with h5py.File(fname, 'w') as f:
_writeMetadata(f)
- _writeLinearTransform(f, xform)
+ _writeAffine(f, xform)
from_ = f.create_group('/From')
to = f.create_group('/To')
- _writeLinearMapping(from_, src)
- _writeLinearMapping(to, ref)
+ _writeSpace(from_, src)
+ _writeSpace(to, ref)
def readNonLinearX5(fname):
@@ -327,8 +549,8 @@ def readNonLinearX5(fname):
with h5py.File(fname, 'r') as f:
field = _readNonLinearTransform(f['/'])
- src = _readLinearMapping(f['/From'])
- ref = _readLinearMapping(f['/To'])
+ src = _readSpace(f['/From'])
+ ref = _readSpace(f['/To'])
# TODO coefficient fields
return nonlinear.DisplacementField(field,
@@ -340,42 +562,7 @@ def readNonLinearX5(fname):
def writeNonLinearX5(fname, field):
"""
- ::
- /Format # "X5"
- /Version # "0.0.1"
- /Metadata # json string containing unstructured metadata
-
- /Transform # the displacement/coefficient field itself
- /Type # "nonlinear"
- /SubType # "displacement" / "deformation"
- /Representation # "cubic bspline" / "quadratic bspline"
- /Inverse # optional pre-calculated inverse
-
- /Pre/Type # "linear"
- /Pre/Transform # ref world-to-[somewhere], to prepare ref
- # world coordinates as inputs to the nonlinear
- # transform
- /Pre/Inverse # optional pre-calculated inverse
- /Post/Type # "linear"
- /Post/Transform # source [somewhere]-to-world, to transform
- # source coordinates produced by the nonlinear
- # transform into world coordinates
- /Post/Inverse # optional pre-calculated inverse
-
- /From/Type # "image"
- /From/Size # voxel dimensions
- /From/Scales # voxel pixdims
- /From/Mapping/Transform # source voxel-to-world sform
- /From/Mapping/Type # "linear"
- /From/Mapping/Inverse # optional inverse
-
- /To/Type # "image"
- /To/Size # voxel dimensions
- /To/Scales # voxel pixdims
- /To/Mapping/Transform # reference voxel-to-world sform
- /To/Mapping/Type # "linear"
- /To/Mapping/Inverse # optional inverse
- """ # noqa
+ """
# TODO coefficient fields