fsl/transform/__init__.py

+#!/usr/bin/env python
+#
+# __init__.py - Functions for working with linear and non-linear FSL
+#               transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with linear and non-linear FSL
+transformations. Functionality is split across the following modules:
+
+.. autosummary::
+   :nosignatures:
+
+   ~fsl.transform.affine
+   ~fsl.transform.flirt
+   ~fsl.transform.fnirt
+   ~fsl.transform.nonlinear
+   ~fsl.transform.x5
+"""

fsl/utils/transform.py → fsl/transform/affine.py

 #!/usr/bin/env python
 #
-# transform.py - Functions for working with affine transformation matrices.
+# affine.py - Utility functions for working with affine transformations.
 #
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
-"""This module provides functions related to 3D image transformations and
-spaces. The following functions are provided:
+"""This module contains utility functions for working with affine
+transformations. The following functions are available:
 
 .. autosummary::
    :nosignatures:
 
    transform
-   transformNormal
    scaleOffsetXform
    invert
+   flip
    concat
    compose
    decompose
@@ -21,9 +21,9 @@ spaces. The following functions are provided:
    rotMatToAxisAngles
    axisAnglesToRotMat
    axisBounds
-   flirtMatrixToSform
-   sformToFlirtMatrix
+   mergeBounds
    rmsdev
+   rescale
 
 And a few more functions are provided for working with vectors:
 
@@ -32,11 +32,13 @@ And a few more functions are provided for working with vectors:
 
    veclength
    normalise
+   transformNormal
 """
 
+
+import collections.abc as abc
 import numpy           as np
 import numpy.linalg    as linalg
-import collections.abc as abc
 
 
 def invert(x):
@@ -60,7 +62,7 @@ def veclength(vec):
 
     Multiple vectors may be passed in, with a shape of ``(n, 3)``.
     """
-    vec = np.array(vec, copy=False).reshape(-1, 3)
+    vec = np.asarray(vec).reshape(-1, 3)
     return np.sqrt(np.einsum('ij,ij->i', vec, vec))
 
 
@@ -69,7 +71,7 @@ def normalise(vec):
 
     Multiple vectors may be passed in, with a shape of ``(n, 3)``.
     """
-    vec = np.array(vec, copy=False).reshape(-1, 3)
+    vec = np.asarray(vec).reshape(-1, 3)
     n   = (vec.T / veclength(vec)).T
 
     if n.size == 3:
@@ -78,7 +80,41 @@ def normalise(vec):
     return n
 
 
-def scaleOffsetXform(scales, offsets):
+def flip(shape, xform, *axes):
+    """Applies a flip/inversion to an affine transform along specified axes.
+
+    :arg shape: The ``(x, y, z)`` shape of the data.
+    :arg xform: Transformation matrix which transforms voxel coordinates
+                to a world coordinate system.
+    :arg axes:  Indices of axes to flip
+    """
+
+    if len(axes) == 0:
+        return xform
+
+    axes = sorted(set(axes))
+    if any(a not in (0, 1, 2) for a in axes):
+        raise ValueError(f'Invalid axis: {axes}')
+
+    los, his = axisBounds(shape, xform, axes, boundary=None)
+    scales   = [1, 1, 1]
+    pre      = [0, 0, 0]
+    post     = [0, 0, 0]
+
+    for ax, lo, hi in zip(axes, los, his):
+        mid        = lo + (hi - lo) / 2
+        scales[ax] = -1
+        pre[   ax] = -mid
+        post[  ax] =  mid
+
+    pre    = scaleOffsetXform(None, pre)
+    post   = scaleOffsetXform(None, post)
+    scales = scaleOffsetXform(scales)
+
+    return concat(post, scales, pre, xform)
+
+
+def scaleOffsetXform(scales=None, offsets=None):
     """Creates and returns an affine transformation matrix which encodes
     the specified scale(s) and offset(s).
 
@@ -94,6 +130,9 @@ def scaleOffsetXform(scales, offsets):
     :returns:     A ``numpy.float32`` array of size :math:`4 \\times 4`.
     """
 
+    if scales  is None: scales  = [1, 1, 1]
+    if offsets is None: offsets = [0, 0, 0]
+
     oktypes = (abc.Sequence, np.ndarray)
 
     if not isinstance(scales,  oktypes): scales  = [scales]
@@ -120,19 +159,28 @@ def scaleOffsetXform(scales, offsets):
     return xform
 
 
-def compose(scales, offsets, rotations, origin=None):
+def compose(scales, offsets, rotations, origin=None, shears=None,
+            scaleAtOrigin=False):
     """Compose a transformation matrix out of the given scales, offsets
     and axis rotations.
 
-    :arg scales:    Sequence of three scale values.
+    :arg scales:        Sequence of three scale values.
+
+    :arg offsets:       Sequence of three offset values.
 
-    :arg offsets:   Sequence of three offset values.
+    :arg rotations:     Sequence of three rotation values, in radians, or
+                        a rotation matrix of shape ``(3, 3)``.
 
-    :arg rotations: Sequence of three rotation values, in radians, or
-                    a rotation matrix of shape ``(3, 3)``.
+    :arg origin:        Origin of rotation.  If not provided, the rotation
+                        origin is ``(0, 0, 0)``.
 
-    :arg origin:    Origin of rotation - must be scaled by the ``scales``.
-                    If not provided, the rotation origin is ``(0, 0, 0)``.
+    :arg shears:        Sequence of three shear values
+
+    :arg scaleAtOrigin: If ``True``, the scaling parameters are applied with
+                        respect to the ``origin``, i.e. so that the location
+                        of ``origin`` is unchanged after scaling.  If
+                        ``False`` (the default), the scaling parameters are
+                        applied with respect to location ``(0, 0, 0)``.
     """
 
     preRotate  = np.eye(4)
@@ -154,6 +202,7 @@ def compose(scales, offsets, rotations, origin=None):
     scale  = np.eye(4, dtype=np.float64)
     offset = np.eye(4, dtype=np.float64)
     rotate = np.eye(4, dtype=np.float64)
+    shear  = np.eye(4, dtype=np.float64)
 
     scale[  0,  0] = scales[ 0]
     scale[  1,  1] = scales[ 1]
@@ -164,10 +213,18 @@ def compose(scales, offsets, rotations, origin=None):
 
     rotate[:3, :3] = rotations
 
-    return concat(offset, postRotate, rotate, preRotate, scale)
+    if shears is not None:
+        shear[0, 1] = shears[0]
+        shear[0, 2] = shears[1]
+        shear[1, 2] = shears[2]
+
+    if scaleAtOrigin:
+        return concat(offset, postRotate, rotate, scale, preRotate, shear)
+    else:
+        return concat(offset, postRotate, rotate, preRotate, scale, shear)
 
 
-def decompose(xform, angles=True):
+def decompose(xform, angles=True, shears=False):
     """Decomposes the given transformation matrix into separate offsets,
     scales, and rotations, according to the algorithm described in:
 
@@ -175,21 +232,24 @@ def decompose(xform, angles=True):
     320-323 in *Graphics Gems II*, James Arvo (editor), Academic Press, 1991,
     ISBN: 0120644819.
 
-    It is assumed that the given transform has no perspective components. Any
-    shears in the affine are discarded.
+    It is assumed that the given transform has no perspective components.
 
-    :arg xform:  A ``(4, 4)`` affine transformation matrix.
+    :arg xform:  A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
 
     :arg angles: If ``True`` (the default), the rotations are returned
                  as axis-angles, in radians. Otherwise, the rotation matrix
                  is returned.
 
+    :arg shears: Defaults to ``False``. If ``True``, shears are returned.
+
     :returns: The following:
 
                - A sequence of three scales
-               - A sequence of three translations
+               - A sequence of three translations (all ``0`` if ``xform``
+                 was a ``(3, 3)`` matrix)
                - A sequence of three rotations, in radians. Or, if
                  ``angles is False``, a rotation matrix.
+               - If ``shears is True``, a sequence of three shears.
     """
 
     # The inline comments in the code below are taken verbatim from
@@ -198,9 +258,13 @@ def decompose(xform, angles=True):
     # The next step is to extract the translations. This is trivial;
     # we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
     # point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
-    xform        = xform.T
-    translations = xform[ 3, :3]
-    xform        = xform[:3, :3]
+    xform = np.array(xform).T
+
+    if xform.shape == (4, 4):
+        translations = xform[ 3, :3]
+        xform        = xform[:3, :3]
+    else:
+        translations = np.array([0, 0, 0])
 
     M1 = xform[0]
     M2 = xform[1]
@@ -209,6 +273,7 @@ def decompose(xform, angles=True):
     # The process of finding the scaling factors and shear parameters
     # is interleaved. First, find s_x = |M'_1|.
     sx = np.sqrt(np.dot(M1, M1))
+    M1 = M1 / sx
 
     # Then, compute an initial value for the xy shear factor,
     # s_xy = M'_1 * M'_2. (this is too large by the y scaling factor).
@@ -225,7 +290,7 @@ def decompose(xform, angles=True):
     # The second row is normalized, and s_xy is divided by s_y to
     # get its final value.
     M2  = M2  / sy
-    sxy = sxy / sy
+    sxy = sxy / sx
 
     # The xz and yz shear factors are computed as in the preceding,
     sxz = np.dot(M1, M3)
@@ -240,8 +305,8 @@ def decompose(xform, angles=True):
     # the third row is normalized, and the xz and yz shear factors are
     # rescaled.
     M3  = M3  / sz
-    sxz = sxz / sz
-    syz = sxz / sz
+    sxz = sxz / sx
+    syz = syz / sy
 
     # The resulting matrix now is a pure rotation matrix, except that it
     # might still include a scale factor of -1. If the determinant of the
@@ -261,7 +326,13 @@ def decompose(xform, angles=True):
     if angles: rotations = rotMatToAxisAngles(R.T)
     else:      rotations = R.T
 
-    return [sx, sy, sz], translations, rotations
+    retval = [np.array([sx, sy, sz]), translations, rotations]
+
+    if shears:
+        retval.append(np.array((sxy, sxz, syz)))
+
+    return tuple(retval)
+
 
 
 def rotMatToAffine(rotmat, origin=None):
@@ -446,11 +517,45 @@ def axisBounds(shape,
     else:      return (lo,    hi)
 
 
+def mergeBounds(*bounds):
+    """Merge multiple ``(lo, hi)`` bounds together. Returns the union of
+    all bounds, i.e. a new set of bounds which encompasses all given bounds.
+
+    Each set of bounds can be provided in one of the following formats:
+
+      - ``((lo, lo, lo), (hi, hi, hi))``
+      - ``((lo, hi), (lo, hi), (lo, hi))``
+
+    The return value will be a tuple of tuples, otherwise having the same
+    format as the inputs.
+    """
+
+    nbounds = len(bounds)
+    bounds  = np.array(bounds)
+
+    if bounds.shape not in ((nbounds, 3, 2), (nbounds, 2, 3)):
+        raise ValueError('Unsupported bounds format')
+
+    if bounds.shape == (nbounds, 2, 3):
+        return ((bounds[:, 0, 0].min(),
+                 bounds[:, 0, 1].min(),
+                 bounds[:, 0, 2].min()),
+                (bounds[:, 1, 0].max(),
+                 bounds[:, 1, 1].max(),
+                 bounds[:, 1, 2].max()))
+    else:
+        return ((bounds[:, 0, 0].min(),
+                 bounds[:, 0, 1].max()),
+                (bounds[:, 1, 0].min(),
+                 bounds[:, 1, 1].max()),
+                (bounds[:, 2, 0].min(),
+                 bounds[:, 2, 1].max()))
+
+
 def transform(p, xform, axes=None, vector=False):
     """Transforms the given set of points ``p`` according to the given affine
     transformation ``xform``.
 
-
     :arg p:      A sequence or array of points of shape :math:`N \\times  3`.
 
     :arg xform:  A ``(4, 4)`` affine transformation matrix with which to
@@ -534,68 +639,6 @@ def _fillPoints(p, axes):
     return newp
 
 
-def flirtMatrixToSform(flirtMat, srcImage, refImage):
-    """Converts the given ``FLIRT`` transformation matrix into a
-    transformation from the source image voxel coordinate system to
-    the reference image world coordinate system.
-
-    FLIRT transformation matrices transform from the source image scaled voxel
-    coordinate system into the reference image scaled voxel coordinate system
-    (voxels scaled by pixdims, with a left-right flip if the image sform has a
-    positive determinant).
-
-    So to construct a transformation from source image voxel coordinates
-    into reference image world coordinates, we need to combine the following:
-
-      1. Source voxels -> Source scaled voxels
-      2. Source scaled voxels -> Reference scaled voxels (the FLIRT matrix)
-      3. Reference scaled voxels -> Reference voxels
-      4. Reference voxels -> Reference world (the reference image sform)
-
-    :arg flirtMat: A ``(4, 4)`` transformation matrix
-    :arg srcImage: Source :class:`.Image`
-    :arg refImage: Reference :class:`.Image`
-    """
-
-    srcScaledVoxelMat    = srcImage.voxToScaledVoxMat
-    refInvScaledVoxelMat = refImage.scaledVoxToVoxMat
-    refVoxToWorldMat     = refImage.voxToWorldMat
-
-    return concat(refVoxToWorldMat,
-                  refInvScaledVoxelMat,
-                  flirtMat,
-                  srcScaledVoxelMat)
-
-
-def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
-    """Under the assumption that the given ``srcImage`` and ``refImage`` share a
-    common world coordinate system (defined by their
-    :attr:`.Nifti.voxToWorldMat` attributes), this function will calculate and
-    return a transformation matrix from the ``srcImage`` scaled voxel
-    coordinate system to the ``refImage`` scaled voxel coordinate system, that
-    can be saved to disk and used with FLIRT, to resample the source image to
-    the reference image.
-
-    :arg srcImage: Source :class:`.Image`
-    :arg refImage: Reference :class:`.Image`
-    :arg srcXform: Optionally used in place of the ``srcImage``
-                   :attr:`.Nifti.voxToWorldMat`
-    """
-
-    srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat
-    srcVoxToWorldMat     = srcImage.voxToWorldMat
-    refWorldToVoxMat     = refImage.worldToVoxMat
-    refVoxToScaledVoxMat = refImage.voxToScaledVoxMat
-
-    if srcXform is not None:
-        srcVoxToWorldMat = srcXform
-
-    return concat(refVoxToScaledVoxMat,
-                  refWorldToVoxMat,
-                  srcVoxToWorldMat,
-                  srcScaledVoxToVoxMat)
-
-
 def rmsdev(T1, T2, R=None, xc=None):
     """Calculates the RMS deviation of the given affine transforms ``T1`` and
     ``T2``. This can be used as a measure of the 'distance' between two
@@ -640,3 +683,59 @@ def rmsdev(T1, T2, R=None, xc=None):
     erms = np.sqrt(erms)
 
     return erms
+
+
+def rescale(oldShape, newShape, origin=None):
+    """Calculates an affine matrix to use for resampling.
+
+    This function generates an affine transformation matrix that can be used
+    to resample an N-D array from ``oldShape`` to ``newShape`` using, for
+    example, ``scipy.ndimage.affine_transform``.
+
+    The matrix will contain scaling factors derived from the ``oldShape /
+    newShape`` ratio, and an offset determined by the ``origin``.
+
+    The default value for ``origin`` (``'centre'``) causes the corner voxel of
+    the output to have the same centre as the corner voxel of the input. If
+    the origin is ``'corner'``, we apply an offset which effectively causes
+    the voxel grid corners of the input and output to be aligned.
+
+    :arg oldShape: Shape of input data
+    :arg newShape: Shape to resample data to
+    :arg origin:   Voxel grid alignment - either ``'centre'`` (the default) or
+                   ``'corner'``
+    :returns:      An affine resampling matrix
+    """
+
+    if origin is None:
+        origin = 'centre'
+
+    oldShape = np.array(oldShape, dtype=float)
+    newShape = np.array(newShape, dtype=float)
+    ndim     = len(oldShape)
+
+    if len(oldShape) != len(newShape):
+        raise ValueError('Shape mismatch')
+
+    # shapes are the same - no rescaling needed
+    if np.all(np.isclose(oldShape, newShape)):
+        return np.eye(ndim + 1)
+
+    # Otherwise we calculate a scaling
+    # matrix from the old/new shape
+    # ratio, and specify an offset
+    # according to the origin
+    ratio = oldShape / newShape
+    scale = np.diag(ratio)
+
+    # Calculate an offset from the origin
+    if   origin == 'centre': offset = [0] * ndim
+    elif origin == 'corner': offset = (ratio - 1) / 2
+
+    # combine the scales and translations
+    # to form thte final affine
+    xform               = np.eye(ndim + 1)
+    xform[:ndim, :ndim] = scale
+    xform[:ndim, -1]    = offset
+
+    return xform

fsl/transform/flirt.py

+#!/usr/bin/env python
+#
+# flirt.py - Functions for working with FLIRT matrices.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with FLIRT affine transformation
+matrices. The following functions are available:
+
+.. autosummary::
+   :nosignatures:
+
+   readFlirt
+   writeFlirt
+   fromFlirt
+   toFlirt
+   flirtMatrixToSform
+   sformToFlirtMatrix
+
+
+
+FLIRT transformation matrices are affine matrices of shape ``(4, 4)`` which
+encode a linear transformation from a source image to a reference image. FLIRT
+matrices are defined in terms of *FSL coordinates*, which is a coordinate
+system where voxels are scaled by pixdims, and with a left-right flip if the
+image ``sform`` has a positive determinant.
+
+
+FLIRT matrices thus encode a transformation from source image FSL coordinates
+to reference image FSL coordinates.
+"""
+
+
+import numpy as np
+
+from .affine import concat
+
+
+def readFlirt(fname):
+    """Reads a FLIRT matrix from a file. """
+    return np.loadtxt(fname)
+
+
+def writeFlirt(xform, fname):
+    """Writes the given FLIRT matrix to a file. """
+    np.savetxt(fname, xform, fmt='%1.15g')
+
+
+def fromFlirt(xform, src, ref, from_='voxel', to='world'):
+    """Convert a FLIRT affine matrix into another affine.
+
+    Given a FLIRT matrix, and the source and reference images with which the
+    FLIRT matrix is associated, converts the matrix into an affine matrix
+    which can transform coordinates from the source image ``from_`` coordinate
+    system to the reference image ``to`` coordinate system.
+
+    The ``from_`` and ``to`` arguments specify the desired spaces that the
+    returned affine should transform between. The default values
+    (``from_='voxel'`` and ``to='world'`` will generate an affine which
+    transforms from voxels in the source image to world-coordinates in the
+    reference image.
+
+    Valid values for the ``from_`` and ``to`` arguments are:
+
+     - ``voxel``: The voxel coordinate system
+     - ``fsl``: The FSL coordinate system (voxels scaled by pixdims, with the
+       X axis inverted if the image sform/qform has a positive determinant)
+     - ``world``  The world coordinate system
+
+    See the :class:`.Nifti` class documentation and the
+    :meth:`.Nifti.getAffine` method for more details.
+
+    :arg xform: ``numpy`` array of shape ``(4, 4)`` containing a FLIRT
+                transformation matrix.
+    :arg src:   :class:`.Nifti` object, the ``xform`` source image
+    :arg ref:   :class:`.Nifti` object, the ``xform`` reference image
+    :arg from_: Desired source coordinate system
+    :arg to:    Desired reference coordinate system
+    :returns:   ``numpy`` array of shape ``(4, 4)`` containing a matrix
+                encoding a transformation from the source ``from_`` to
+                the reference ``to`` coordinate systems.
+    """
+    premat  = src.getAffine(from_, 'fsl')
+    postmat = ref.getAffine('fsl', to)
+    return concat(postmat, xform, premat)
+
+
+def toFlirt(xform, src, ref, from_='voxel', to='world'):
+    """Convert an affine matrix into a FLIRT matrix.
+
+    :returns:   ``numpy`` array of shape ``(4, 4)`` containing a matrix
+                encoding a transformation from the source ``from_`` to
+                the reference ``to`` coordinate systems.
+    :arg src:   :class:`.Nifti` object, the ``xform`` source image
+    :arg ref:   :class:`.Nifti` object, the ``xform`` reference image
+    :arg from_: ``xform`` source coordinate system
+    :arg to:    ``xform`` target coordinate system
+    :returns:   A ``numpy`` array of shape ``(4, 4)`` containing a FLIRT
+                transformation matrix from ``src`` to ``ref``.
+    """
+    premat  = src.getAffine('fsl', from_)
+    postmat = ref.getAffine(to, 'fsl')
+    return concat(postmat, xform, premat)
+
+
+def flirtMatrixToSform(flirtMat, srcImage, refImage):
+    """Converts the given ``FLIRT`` transformation matrix into a
+    transformation from the source image voxel coordinate system to
+    the reference image world coordinate system.
+
+    To construct a transformation from source image voxel coordinates into
+    reference image world coordinates, we need to combine the following:
+
+      1. Source voxels -> Source scaled voxels
+      2. Source scaled voxels -> Reference scaled voxels (the FLIRT matrix)
+      3. Reference scaled voxels -> Reference voxels
+      4. Reference voxels -> Reference world (the reference image sform)
+
+    :arg flirtMat: A ``(4, 4)`` transformation matrix
+    :arg srcImage: Source :class:`.Image`
+    :arg refImage: Reference :class:`.Image`
+    """
+    return fromFlirt(flirtMat, srcImage, refImage, 'voxel', 'world')
+
+
+def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
+    """Under the assumption that the given ``srcImage`` and ``refImage`` share a
+    common world coordinate system (defined by their
+    :attr:`.Nifti.voxToWorldMat` attributes), this function will calculate and
+    return a transformation matrix from the ``srcImage`` FSL coordinate system
+    to the ``refImage`` FSL coordinate system, that can be saved to disk and
+    used with FLIRT, to resample the source image to the space of the
+    reference image.
+
+    :arg srcImage: Source :class:`.Image`
+    :arg refImage: Reference :class:`.Image`
+    :arg srcXform: Optionally used in place of the ``srcImage``
+                   :attr:`.Nifti.voxToWorldMat`
+    """
+
+    srcFSLToVoxMat   = srcImage.getAffine('fsl',   'voxel')
+    srcVoxToWorldMat = srcImage.getAffine('voxel', 'world')
+    refWorldToVoxMat = refImage.getAffine('world', 'voxel')
+    refVoxToFSLMat   = refImage.getAffine('voxel', 'fsl')
+
+    if srcXform is not None:
+        srcVoxToWorldMat = srcXform
+
+    return concat(refVoxToFSLMat,
+                  refWorldToVoxMat,
+                  srcVoxToWorldMat,
+                  srcFSLToVoxMat)

fsl/transform/fnirt.py

+#!/usr/bin/env python
+#
+# fnirt.py - Functions for working with FNIRT non-linear transformations.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module contains functions for working with FNIRT non-linear
+transformations. The following functions are available:
+
+.. autosummary::
+   :nosignatures:
+
+   readFnirt
+   toFnirt
+   fromFnirt
+
+
+Non-linear registration using FNIRT
+-----------------------------------
+
+
+FNIRT is used to calculate a non-linear registration from a source image to a
+reference image. FNIRT outputs the resulting non-linear transformation as
+either:
+
+ - A deformation/warp field which contains displacements or coordinates.
+ - A coefficient field which can be used to generate a warp field.
+
+
+Non-linear registration using FNIRT generally follows the process depicted
+here:
+
+
+.. image:: images/nonlinear_registration_process.png
+   :width: 80%
+   :align: center
+
+
+First, an initial linear registration is performed from the source image to
+the reference image using FLIRT; this provides an initial global alignment
+which can be used as the starting point for the non-linear registration. Next,
+FNIRT is used to non-linearly register the aligned source image to the
+reference image.  Importantly, both of these steps are performed using FSL
+coordinates.
+
+
+Note that we have three spaces, and three sets of coordinate systems, to
+consider:
+
+ 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
+
+ 3. Reference image space
+
+
+The initial affine registration calculates a transformation between spaces 1
+and 2, and the non-linear registration calculates a transformation between
+spaces 2 and 3. Note that the fields-of-view for spaces 2 and 3 are
+equivalent.
+
+
+The non-linear transformation file generated by FNIRT will contain the initial
+linear registration, with it either being encoded directly into the warps (for
+a warp field), or being stored in the NIfTI header (for a coefficient field).
+
+
+FNIRT warp fields
+^^^^^^^^^^^^^^^^^
+
+
+A FNIRT deformation field (a.k.a. warp field) is defined in the same space as
+the reference image, and may contain:
+
+ - *relative displacements*, where each voxel in the warp 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 coordinates*, where each voxel in the warp field simply contains
+   the source image coordinates which correspond to those reference image
+   coordinates.
+
+.. note:: FNIRT deformation field files give no indication as to whether they
+          contain relative displacements or absolute coordinates, so heuristics
+          must be used to infer what is stored in a a particular file. The
+          :func:`.nonlinear.detectDeformationType` function can be used to
+          determine whether a file contains relative displacements or absolute
+          coordinates.
+
+If an initial linear registration was used as the starting point for FNIRT,
+this is encoded into the displacements/coordinates themselves, so they can be
+used to transform from the reference image to the *original* source image
+space.
+
+
+.. image:: images/fnirt_warp_field.png
+   :width: 80%
+   :align: center
+
+
+FNIRT coefficient fields
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+
+A FNIRT coefficient field 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 added to that location's
+reference image coordinates, in order to determine the corresponding source
+image coordinates.
+
+
+If an initial linear registration was used as the starting point for FNIRT,
+the generated displacement field will encode a transformation to *aligned*
+source image coordinates, and the initial affine will be stored in the NIfTI
+header of the coefficient field file.
+
+
+.. image:: images/fnirt_coefficient_field.png
+   :width: 80%
+   :align: center
+"""
+
+
+import logging
+
+import nibabel as nib
+import numpy   as np
+
+import fsl.data.constants as constants
+import fsl.data.image     as fslimage
+from . import                affine
+from . import                nonlinear
+
+
+log = logging.getLogger(__name__)
+
+
+def _readFnirtDeformationField(fname, img, src, ref, defType=None):
+    """Loads ``fname``, assumed to be a FNIRT deformation field.
+
+    :arg fname:   File name of FNIRT deformation field
+
+    :arg img:     ``fname`` loaded as an :class:`.Image`
+
+    :arg src:     Source image
+
+    :arg ref:     Reference image
+
+    :arg defType: Deformation type - either ``'absolute'`` or ``'relative'``.
+                  If not provided, is automatically inferred from the data.
+
+    :return:      A :class:`.DeformationField` object
+    """
+    return nonlinear.DeformationField(fname,
+                                      src,
+                                      ref,
+                                      srcSpace='fsl',
+                                      refSpace='fsl',
+                                      defType=defType)
+
+
+def _readFnirtCoefficientField(fname, img, src, ref):
+    """Loads ``fname``, assumed to be a FNIRT coefficient field.
+
+    :arg fname: File name of FNIRT coefficient field
+    :arg img:   ``fname`` loaded as an :class:`.Image`
+    :arg src:   Source image
+    :arg ref:   Reference image
+    :return:    A :class:`.CoefficientField`
+    """
+
+    # FNIRT uses NIFTI header fields in
+    # non-standard ways to store some
+    # additional information about the
+    # coefficient field. See
+    # $FSLDIR/src/fnirt/fnirt_file_writer.cpp
+    # for more details.
+
+    # The field type (quadratic, cubic,
+    # or discrete-cosine-transform) is
+    # inferred from the intent. There is
+    # no support in this implementation
+    # for DCT fields
+    cubics = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS)
+    quads  = (constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
+
+    if img.intent in cubics:
+        fieldType = 'cubic'
+    elif img.intent in quads:
+        fieldType = 'quadratic'
+    else:
+        fieldType = 'cubic'
+        log.warning('Unrecognised/unsupported coefficient '
+                    'field type (assuming cubic b-spline): '
+                    '{}'.format(img.intent))
+
+    # Knot spacing (in voxels) is
+    # stored in the pixdims
+    knotSpacing = img.pixdim[:3]
+
+    # The sform contains an initial
+    # global src-to-ref affine
+    # (the starting point for the
+    # non-linear registration). This
+    # is encoded as a flirt matrix,
+    # i.e. it transforms from
+    # source-scaled-voxels to
+    # ref-scaled-voxels
+    srcToRefMat = img.header.get_sform()
+
+    # The fieldToRefMat affine tells
+    # the CoefficientField class how
+    # to transform coefficient field
+    # voxel coordinates into
+    # displacement field/reference
+    # image voxel coordinates.
+    fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
+
+    # But if the provided reference has
+    # different resolution to the
+    # reference that was originally
+    # used to generate the warp field,
+    # we need to adjust the field
+    # accordingly. We assume that the
+    # references are aligned in the FSL
+    # coordinate system, so simply apply
+    # a scaling factor calculated by
+    # dividing the original reference
+    # pixdims by the provided reference
+    # pixdims.
+    refPixdims = np.array([img.header['intent_p1'],
+                           img.header['intent_p2'],
+                           img.header['intent_p3']])
+
+    if not np.all(np.isclose(ref.pixdim[:3], refPixdims)):
+        fieldToRefMat = affine.concat(
+            affine.scaleOffsetXform(refPixdims / ref.pixdim[:3], 0),
+            fieldToRefMat)
+
+    return nonlinear.CoefficientField(fname,
+                                      src,
+                                      ref,
+                                      srcSpace='fsl',
+                                      refSpace='fsl',
+                                      fieldType=fieldType,
+                                      knotSpacing=knotSpacing,
+                                      srcToRefMat=srcToRefMat,
+                                      fieldToRefMat=fieldToRefMat)
+
+
+def readFnirt(fname, src, ref, defType=None, intent=None):
+    """Reads a non-linear FNIRT transformation image, returning
+    a :class:`.DeformationField` or :class:`.CoefficientField` depending
+    on the file type.
+
+    :arg fname:    File name of FNIRT transformation
+    :arg src:      Source image
+    :arg ref:      Reference image
+    :arg defType:  Deformation type - either ``'absolute'`` or ``'relative'``.
+                   Only used if the file is a deformation field. If not
+                   provided, is automatically inferred from the data.
+    :arg intent:   NIFTI intent code of ``fname``. e.g.
+                   :attr:`.constants.FSL_FNIRT_DISPLACEMENT_FIELD`. If not
+                   provided, the intent is read from the image header.
+    """
+
+    if defType not in (None, 'absolute', 'relative'):
+        raise ValueError('defType must be None, "absolute" or "relative" '
+                         '(passed in as {})'.format(defType))
+
+    # Figure out whether the file is a
+    # deformation field or a coefficient
+    # field by checking the intent code.
+    # If the intent is provided, assume
+    # that the caller knows the type of
+    # the field.
+    img    = fslimage.Image(fname)
+    intent = intent or img.intent
+    disps  = (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
+              constants.FSL_TOPUP_FIELD)
+    coefs  = (constants.FSL_CUBIC_SPLINE_COEFFICIENTS,
+              constants.FSL_DCT_COEFFICIENTS,
+              constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_CUBIC_SPLINE_COEFFICIENTS,
+              constants.FSL_TOPUP_QUADRATIC_SPLINE_COEFFICIENTS)
+
+    if intent in disps:
+        return _readFnirtDeformationField(fname, img, src, ref, defType)
+    elif intent in coefs:
+        return _readFnirtCoefficientField(fname, img, src, ref)
+    else:
+        raise ValueError('Cannot determine type of nonlinear warp field '
+                         '{} (intent code: {})'.format(fname, intent))
+
+
+def toFnirt(field):
+    """Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
+    :class:`.DeformationField` or :class:`.CoefficientField`.
+
+    :arg field: :class:`.NonLinearTransform` to convert
+    :return:    A FNIRT-compatible :class:`.DeformationField` or
+                :class:`.CoefficientField`.
+    """
+
+    # If we have a coefficient field
+    # which transforms between fsl
+    # space, we can just create a copy.
+    if isinstance(field, nonlinear.CoefficientField) and \
+       (field.srcSpace == 'fsl' and field.refSpace == 'fsl'):
+
+        # We start with a nibabel image,
+        # as we need to mess with the header
+        # fields to store all of the FNIRT
+        # coefficient field information
+        fieldBase = nib.nifti1.Nifti1Image(field.data, None)
+
+        # Set the intent
+        if field.fieldType == 'cubic':
+            intent = constants.FSL_CUBIC_SPLINE_COEFFICIENTS
+        elif field.fieldType == 'quadratic':
+            intent = constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS
+        fieldBase.header['intent_code'] = intent
+
+        # Reference image pixdims are
+        # stored in the intent code
+        # parameters.
+        fieldBase.header['intent_p1'] = field.ref.pixdim[0]
+        fieldBase.header['intent_p2'] = field.ref.pixdim[1]
+        fieldBase.header['intent_p3'] = field.ref.pixdim[2]
+
+        # Pixdims are used to
+        # store the knot spacing,
+        pixdims      = list(field.knotSpacing) + [1]
+        qform        = np.diag(pixdims)
+
+        # The sform is used to store the
+        # initial src-to-ref affine
+        if field.srcToRefMat is not None: sform = field.srcToRefMat
+        else:                             sform = np.eye(4)
+
+        # The qform offsets are
+        # used to store the
+        # reference image shape
+        qform[:3, 3] = field.ref.shape[:3]
+
+        fieldBase.header.set_zooms(pixdims)
+        fieldBase.set_sform(sform, 1)
+        fieldBase.set_qform(qform, 1)
+        fieldBase.update_header()
+
+        field = nonlinear.CoefficientField(
+            fieldBase,
+            src=field.src,
+            ref=field.ref,
+            srcSpace='fsl',
+            refSpace='fsl',
+            fieldType=field.fieldType,
+            knotSpacing=field.knotSpacing,
+            fieldToRefMat=field.fieldToRefMat,
+            srcToRefMat=field.srcToRefMat)
+
+    # Otherwise we have a non-FSL coefficient
+    # field, or a deformation field.
+    else:
+
+        # We can't convert a CoefficientField
+        # which doesn't transform in FSL
+        # coordinates, because the coefficients
+        # will have been calculated between some
+        # other source/reference coordinate
+        # systems, and we can't adjust the
+        # coefficients to encode an FSL->FSL
+        # deformation.
+        if isinstance(field, nonlinear.CoefficientField):
+            field = nonlinear.coefficientFieldToDeformationField(field)
+
+        # Again, if we have a displacement
+        # field which transforms between
+        # fsl spaces, we can just take a copy
+        if field.srcSpace == 'fsl' and field.refSpace == 'fsl':
+            field = nonlinear.DeformationField(
+                field.data,
+                header=field.header,
+                src=field.src,
+                ref=field.ref,
+                defType=field.deformationType)
+
+        # Otherwise we have to adjust the
+        # displacements so they transform
+        # between fsl coordinates.
+        field = nonlinear.convertDeformationSpace(
+            field, from_='fsl', to='fsl')
+
+        field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
+
+    return field
+
+
+def fromFnirt(field, from_='world', to='world'):
+    """Convert a FNIRT-style :class:`.NonLinearTransform` to a generic
+    :class:`.DeformationField`.
+
+    :arg field: A FNIRT-style :class:`.CoefficientField` or
+                :class:`.DeformationField`
+
+    :arg from_: Desired reference image coordinate system
+
+    :arg to:    Desired source image coordinate system
+
+    :return:    A :class:`.DeformationField` which contains displacements
+                from the reference image ``from_`` cordinate system to the
+                source image ``to`` coordinate syste.
+    """
+    if isinstance(field, nonlinear.CoefficientField):
+        field = nonlinear.coefficientFieldToDeformationField(field)
+    return nonlinear.convertDeformationSpace(field, from_=from_, to=to)