DOC: Notes on X5 formats. Not finished. Other minor doc additions

fsl/scripts/fsl_convert_x5.py

 #!/usr/bin/env python
 #
-# fsl_convert_x5.py -
+# fsl_convert_x5.py - Convert between FSL and X5 transformation files.
 #
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
+"""This script can be used to convert between FSL and X5 linear and non-linear
+transformation file formats.
+"""
 
 
 import os.path as op
@@ -78,12 +81,16 @@ def flirtToX5(args):
 
 
 def X5ToFlirt(args):
+    """Convert a linear X5 transformation file to a FLIRT matrix. """
     xform, src, ref = transform.readLinearX5(args.input)
     xform           = transform.toFlirt(xform, src, ref, 'world', 'world')
     transform.writeFlirt(xform, args.output)
 
 
 def fnirtToX5(args):
+    """Convert a non-linear FNIRT transformation into an X5 transformation
+    file.
+    """
     src   = fslimage.Image(args.source,    loadData=False)
     ref   = fslimage.Image(args.reference, loadData=False)
     field = transform.readFnirt(args.input,
@@ -95,12 +102,16 @@ def fnirtToX5(args):
 
 
 def X5ToFnirt(args):
+    """Convert a non-linear X5 transformation file to a FNIRT-style
+    transformation file.
+    """
     field = transform.readNonLinearX5(args.input)
     field = transform.toFnirt(field, 'world', 'world')
     transform.writeFnirt(field, args.output)
 
 
 def doFlirt(args):
+    """Converts a linear FIRT transformation file to or from the X5 format. """
     infmt  = args.input_format
     outfmt = args.output_format
 
@@ -110,6 +121,9 @@ def doFlirt(args):
 
 
 def doFnirt(args):
+    """Converts a non-linear FNIRT transformation file to or from the X5
+    format.
+    """
     infmt  = args.input_format
     outfmt = args.output_format
 
@@ -119,15 +133,22 @@ def doFnirt(args):
 
 
 def main(args=None):
+    """Entry point. Calls :func:`doFlirt` or :func:`doFnirt` depending on
+    the sub-command specified in the arguments.
+
+    :arg args: Sequence of command-line arguments. If not provided,
+               ``sys.argv`` is used.
+    """
 
     if args is None:
         args = sys.argv[1:]
 
-    args   = parseArgs(args)
-    ctype  = args.ctype
+    args = parseArgs(args)
+
+    if   args.ctype == 'flirt': doFlirt(args)
+    elif args.ctype == 'fnirt': doFnirt(args)
 
-    if   ctype == 'flirt': doFlirt(args)
-    elif ctype == 'fnirt': doFnirt(args)
+    return 0
 
 
 if __name__ == '__main__':

fsl/transform/x5.py

 #!/usr/bin/env python
 #
-# x5.py -
+# x5.py - Functions for working with BIDS X5 files.
 #
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
-"""
+"""This module contains functions for reading/writing linear/non-linear FSL
+transformations from/to BIDS X5 files.
+
+
+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
+an X5 file, the source image is referred to as the ``From`` space, and the
+reference image the ``To`` space.
+
+
+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.
+
+
+To simplify the file format efinitions 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.
+
+
+*affine*
+^^^^^^^^
+
+
+A HDF5 group which is listed as being of type "affine" contains an affine
+transformation, which can be used to transform coordinates from one space into
+another. Groups of type "affine" have the following fields:
+
+
++---------------+-----------+-------------------------------------------------+
+| Name          | Type      | Value/Description                               |
++---------------+-----------+-------------------------------------------------+
+| ``Type``      | attribute | ``'linear'``                                    |
+| ``Transform`` | dataset   | The transformation itself - a 4x4 ``float64``   |
+|               |           | affine transformation                           |
+| ``Inverse``   | dataset   | Optional pre-calculated inverse                 |
++---------------+-----------+-------------------------------------------------+
+
+
+*mapping*
+^^^^^^^^^
+
+
+A HDF5 group which is listed as being of type "mapping" 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.
+
+
+Groups of type "mapping" have the following fields:
+
+
++-------------+-----------+---------------------------------------------------+
+| 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")                                          |
++-------------+-----------+---------------------------------------------------+
+
+
+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**.
+
+
+Linear X5 transformation files are assumed to adhere to the HDF5 structure
+defined in the table below.  All fields are required.
+
+
++-----------------------------+-----------+-----------------------------------+
+| 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 |
++-----------------------------+-----------+-----------------------------------+
+
+
+Non-linear X5 transformation files
+----------------------------------
+
+
+Non-linear X5 transformation files contain a non-linear transformation between
+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.
+
+ - A quadratic or cubic B-spline *coefficient field*, which contains
+   coefficients defined in a coarse grid which overlays 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)
+
+ - *absolute displacement* which simply contain the source image coordinates
+   at each voxel in the reference image space.
+
+
+Displacement field transformations define displacements between the refernece
+image *world* coordinate system, and the source image *world* coordinate
+system.
+
+
+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.
+
+
+Non-linear X5 transformation files are assumed to adhere to the following
+HDF5 structure::
+
+
++-----------------------------+-----------+-----------------------------------+
+| Name                        | Type      | Value/Description                 |
++-----------------------------+-----------+-----------------------------------+
+| *Metadata* - same as for linear transformation files                        |
++-----------------------------------------------------------------------------+
+| *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           |
+|-----------------------------+-----------+-----------------------------------+
+| *Coefficient field parameters* (required for ``'coefficient'`` files)       |
+|-----------------------------+-----------+-----------------------------------+
+| ``/Parameters/Spacing``          | attribute |
+| ``/Parameters/ReferenceToField`` | affine    |    | attribute |
++-----------------------------+-----------+-----------------------------------+
 
-Functions for reading/writing linear/non-linear FSL transformations from/to
-BIDS X5 files.
 """
 
+
 import json
 
-import numpy        as np
-import numpy.linalg as npla
-import nibabel      as nib
+import numpy   as np
+import nibabel as nib
 import h5py
 
 import fsl.version as version
+from . import         affine
+
+
+X5_FORMAT  = 'X5'
+X5_VERSION = '0.0.1'
 
 
 def _writeMetadata(group):
-    group.attrs['Format']   = 'X5'
-    group.attrs['Version']  = '0.0.1'
+    """Writes a metadata block into the given group.
+
+    :arg group: A ``h5py.Group`` object.
+    """
+    group.attrs['Format']   = X5_FORMAT
+    group.attrs['Version']  = X5_VERSION
     group.attrs['Metadata'] = json.dumps({'fslpy' : version.__version__})
 
 
 def _readLinearTransform(group):
+    """Reads a linear transformation from the given group,
+
+    :arg group: A ``h5py.Group`` object.
+    :returns:   ``numpy`` array containing a ``(4, 4)`` affine transformation.
+    """
+
     if group.attrs['Type'] != 'linear':
         raise ValueError('Not a linear transform')
-    return np.array(group['Transform'])
+
+    xform = np.array(group['Transform'])
+
+    if xform.shape != (4, 4):
+        raise ValueError('Not a linear transform')
+
+    return xform
 
 
 def _writeLinearTransform(group, xform):
+    """Writes the given linear transformation and its inverse to the given
+    group.
+
+    :arg group: A ``h5py.Group`` object.
+    :arg xform: ``numpy`` array containing a ``(4, 4)`` affine transformation.
+    """
 
-    xform = np.asarray(xform,           dtype=np.float32)
-    inv   = np.asarray(npla.inv(xform), dtype=np.float32)
+    xform = np.asarray(xform,                 dtype=np.float64)
+    inv   = np.asarray(affine.inverse(xform), dtype=np.float64)
 
     group.attrs['Type'] = 'linear'
     group.create_dataset('Transform', data=xform)
@@ -43,6 +243,12 @@ def _writeLinearTransform(group, xform):
 
 
 def _readLinearMapping(group):
+    """Reads a linear mapping, defining a source or reference space, from the
+    given group.
+
+    :arg group: A ``h5py.Group`` object.
+    :returns:   :class:`.Nifti` object. defining the mapping
+    """
 
     import fsl.data.image as fslimage
 
@@ -61,6 +267,11 @@ def _readLinearMapping(group):
 
 
 def _writeLinearMapping(group, img):
+    """Writes a linear mapping specified by ``img`` to the given group.
+
+    :arg group: A ``h5py.Group`` object.
+    :arg img:   :class:`.Nifti` object. defining the mapping
+    """
     group.attrs['Type']   = 'image'
     group.attrs['Size']   = np.asarray(img.shape[ :3], np.uint32)
     group.attrs['Scales'] = np.asarray(img.pixdim[:3], np.float32)
@@ -95,33 +306,7 @@ def readLinearX5(fname):
 
 def writeLinearX5(fname, xform, src, ref):
     """
-
-
-    ::
-        /Format                       # "X5"
-        /Version                      # "0.0.1"
-        /Metadata                     # json string containing unstructured metadata
-
-        /Type                         # "linear"
-        /Transform                    # the transform itself
-        /Inverse                      # optional pre-calculated inverse
-
-        /From/Type                    # "image" - could in principle be something other than
-                                      # "image" (e.g. "surface"), in which case the "Size" and
-                                      # "Scales" entries might be replaced with something else
-        /From/Size                    # voxel dimensions
-        /From/Scales                  # voxel pixdims
-        /From/Mapping/Type            # "linear"
-        /From/Mapping/Transform       # source voxel-to-world sform
-        /From/Mapping/Inverse         # optional inverse
-
-        /To/Type                      # "image"
-        /To/Size                      # voxel dimensions
-        /To/Scales                    # voxel pixdims
-        /To/Mapping/Type              # "linear"
-        /To/Mapping/Transform         # reference voxel-to-world sform
-        /To/Mapping/Inverse           # optional inverse
-    """  # noqa
+    """
 
     with h5py.File(fname, 'w') as f:
         _writeMetadata(f)