RF: Move FSL-specific coefficient field stuff from nonlinear.py into fnirt.py

- CoefficientField can now represent a generic b-spline coefficient
field. Cleaned up fnirt.py, other tweaks and comments

fsl/transform/__init__.py

@@ -36,7 +36,6 @@ from .flirt  import (  # noqa
 
 from .fnirt import (  # noqa
     readFnirt,
-    writeFnirt,
     toFnirt,
     fromFnirt)
 

fsl/transform/fnirt.py

@@ -11,7 +11,8 @@ transformation matrices. The following functions are available:
    :nosignatures:
 
    readFnirt
-   writeFnirt
+   toFnirt
+   fromFnirt
 """
 
 
@@ -23,59 +24,167 @@ import fsl.data.constants as constants
 log = logging.getLogger(__name__)
 
 
-def readFnirt(fname, src, ref, dispType=None):
+def _readFnirtDisplacementField(fname, img, src, ref):
+    """Loads ``fname``, assumed to be a FNIRT displacement field.
+
+    :arg fname: File name of FNIRT displacement field
+    :arg img:   ``fname`` loaded as an :class:`.Image`
+    :arg src:   Source image
+    :arg ref:   Reference image
+    :return:    A :class:`.DisplacementField`
     """
+    from . import nonlinear
+    return nonlinear.DisplacementField(fname,
+                                       src,
+                                       ref,
+                                       srcSpace='fsl',
+                                       refSpace='fsl')
+
+
+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`
+    """
+
+    from . import affine
+    from . import nonlinear
+
+    # 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)
+    srcToRefMat = img.header.get_sform()
+
+    # The fieldToRefMat affine allows us
+    # to transform coefficient field voxel
+    # coordinates into displacement field/
+    # reference image voxel coordinates.
+    fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
+
+    return nonlinear.CoefficientField(fname,
+                                      src,
+                                      ref,
+                                      srcSpace='fsl',
+                                      refSpace='fsl',
+                                      fieldType=fieldType,
+                                      knotSpacing=knotSpacing,
+                                      srcToRefMat=srcToRefMat,
+                                      fieldToRefMat=fieldToRefMat)
+
+
+def readFnirt(fname, src, ref):
+    """Reads a non-linear FNIRT transformation image, returning
+    a :class:`.DisplacementField` or :class:`.CoefficientField` depending
+    on the file type.
+
+    :arg fname: File name of FNIRT transformation
+    :arg src:   Source image
+    :arg ref:   Reference image
     """
 
     # Figure out whether the file
     # is a displacement field or
     # a coefficient field
-    import fsl.data.image     as fslimage
-    from   .              import nonlinear
+    import fsl.data.image as fslimage
+
+    img   = fslimage.Image(fname, loadData=False)
+    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 img.intent in disps:
+        return _readFnirtDisplacementField(fname, img, src, ref)
+    elif img.intent in coefs:
+        return _readFnirtCoefficientField(fname, img, src, ref)
+    else:
+        raise ValueError('Cannot determine type of nonlinear '
+                         'file {}'.format(fname))
 
-    img = fslimage.Image(fname, loadData=False)
 
-    dispfields = (constants.FSL_FNIRT_DISPLACEMENT_FIELD,
-                  constants.FSL_TOPUP_FIELD)
-    coeffields = (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)
+def toFnirt(field):
+    """Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
+    :class:`.DisplacementField`.
 
-    kwargs = {
-        'src'      : src,
-        'ref'      : ref,
-        'srcSpace' : 'fsl',
-        'refSpace' : 'fsl',
-        'dispType' : None,
-    }
+    :arg field: :class:`.NonLinearTransform` to convert
+    :return:    A FNIRT-compatible :class:`.DisplacementField`.
+    """
 
-    if img.intent in dispfields:
-        return nonlinear.DisplacementField(fname, **kwargs)
+    from . import nonlinear
 
-    elif img.intent in coeffields:
-        pass  # return nonlinear.CoefficientField(fname, **kwargs)
+    # We can't convert a CoefficientField,
+    # 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.coefficientFieldToDisplacementField(field)
 
-    else:
-        raise ValueError('Cannot determine type of nonlinear '
-                         'file {}'.format(fname))
+    field = nonlinear.convertDisplacementSpace(field, from_='fsl', to='fsl')
+    field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
 
+    return field
 
-def writeFnirt(field, fname):
-    """
-    """
-    field.save(fname)
 
+def fromFnirt(field, from_='world', to='world'):
+    """Convert a FNIRT-style :class:`.NonLinearTransform` to a generic
+    :class:`.DisplacementField`.
 
-def toFnirt(field):
-    pass
+    :arg field: A FNIRT-style :class:`.CoefficientField` or
+                :class:`.DisplacementField`
 
+    :arg from_: Desired reference image coordinate system
 
-def fromFnirt(field, from_='voxel', to='world'):
-    """
-    """
+    :arg to:    Desired source image coordinate system
 
+    :return:    A :class:`.DisplacementField` which contains displacements
+                from the reference image ``from_`` cordinate system to the
+                source image ``to`` coordinate syste.
+    """
     from . import nonlinear
 
+    # see comments in toFnirt
+    if isinstance(field, nonlinear.CoefficientField):
+        field = nonlinear.coefficientFieldToDisplacementField(field)
+
     return nonlinear.convertDisplacementSpace(field, from_=from_, to=to)

fsl/transform/nonlinear.py

@@ -27,8 +27,7 @@ import itertools as it
 
 import numpy as np
 
-import fsl.data.constants as constants
-import fsl.data.image     as fslimage
+import fsl.data.image as fslimage
 
 from . import affine
 
@@ -74,10 +73,10 @@ class NonLinearTransform(fslimage.Image):
         :arg ref:      :class:`.Nifti` representing the reference image.
 
         :arg srcSpace: Coordinate system in the source image that this
-                       ``NonLinearTransform`` maps from. Defaults to ``'fsl'``.
+                       ``NonLinearTransform`` maps to. Defaults to ``'fsl'``.
 
         :arg refSpace: Coordinate system in the reference image that this
-                       ``NonLinearTransform`` maps to. Defaults to ``'fsl'``.
+                       ``NonLinearTransform`` maps from. Defaults to ``'fsl'``.
 
         All other arguments are passed through to :meth:`.Image.__init__`.
         """
@@ -275,53 +274,41 @@ class CoefficientField(NonLinearTransform):
     generated by FNIRT.
     """
 
-    def __init__(self, image, src, ref, **kwargs):
+    def __init__(self,
+                 image,
+                 src,
+                 ref,
+                 srcSpace,
+                 refSpace,
+                 fieldType,
+                 knotSpacing,
+                 srcToRefMat,
+                 fieldToRefMat,
+                 **kwargs):
         """Create a ``CoefficientField``.
 
-        :arg image:
-        :arg src:
-        :arg ref:
+        :arg fieldType:
+        :arg knotSpacing:
+        :arg srcToRefMat:
+        :arg fieldToRefMat:
         """
 
-        NonLinearTransform.__init__(self, image, src, ref, **kwargs)
-
-        # 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
-        if self.intent == constants.FSL_CUBIC_SPLINE_COEFFICIENTS:
-            self.__fieldType = 'cubic'
-        elif self.intent == constants.FSL_QUADRATIC_SPLINE_COEFFICIENTS:
-            self.__fieldType = 'quadratic'
-        else:
-            self.__fieldType = 'cubic'
-            log.warning('Unrecognised/unsupported coefficient field type '
-                        '(assuming cubic): {}'.format(self.intent))
-
-        # Knot spacing (in voxels) is
-        # stored in the pixdims
-        kx, ky, kz = self.pixdim[:3]
-        self.__knotSpacing = (kx, ky, kz)
-
-        # The sform contains an initial
-        # global src-to-ref affine
-        # (the starting point for the
-        # non-linear registration)
-        self.__srcToRefMat = self.header.get_sform()
-
-        # The fieldToRefMat affine allows us
-        # to transform coefficient field voxel
-        # coordinates into displacement field/
-        # reference image voxel coordinates.
-        self.__fieldToRefMat = affine.scaleOffsetXform((kx, ky, kz), 0)
+        if fieldType not in ('quadratic', 'cubic'):
+            raise ValueError('Unsupported field type: {}'.format(fieldType))
+
+        NonLinearTransform.__init__(self,
+                                    image,
+                                    src,
+                                    ref,
+                                    srcSpace,
+                                    refSpace,
+                                    refSpace='voxel',
+                                    **kwargs)
+
+        self.__fieldType     = fieldType
+        self.__knotSpacing   = tuple(knotSpacing)
+        self.__srcToRefMat   = np.copy(srcToRefMat)
+        self.__fieldToRefMat = np.copy(fieldToRefMat)
         self.__refToFieldMat = affine.invert(self.__fieldToRefMat)
 
 
@@ -348,6 +335,7 @@ class CoefficientField(NonLinearTransform):
         """
         return self.__knotSpacing
 
+
     @property
     def fieldToRefMat(self):
         """Return an affine transformation which can transform coefficient
@@ -580,12 +568,13 @@ def convertDisplacementSpace(field, from_, to):
 
 
 def coefficientFieldToDisplacementField(field):
-    """Convert a FNIRT quadratic or cubic B-spline coefficient field into
-    a relative displacement field.
+    """Convert a :class:`CoefficientField` into a relative
+    :class:`DisplacementField`.
 
     :arg field: :class:`CoefficientField` to convert
     :return:    :class:`DisplacementField` calculated from ``field``
     """
+
     ix, iy, iz = field.ref.shape[:3]
     x,  y,  z  = np.meshgrid(np.arange(ix),
                              np.arange(iy),
@@ -596,4 +585,10 @@ def coefficientFieldToDisplacementField(field):
     xyz        = np.vstack((x, y, z)).T
     disps      = field.displacements(xyz).reshape((ix, iy, iz, 3))
 
-    return DisplacementField(disps, field.src, field.ref, dispType='relative')
+    return DisplacementField(disps,
+                             src=field.src,
+                             ref=field.ref,
+                             srcSpace=field.srcSpace,
+                             refSpace=field.refSpace,
+                             header=field.ref.header,
+                             dispType='relative')

fsl/transform/x5.py

@@ -111,7 +111,7 @@ def writeLinearX5(fname, xform, src, ref):
                                       # "Scales" entries might be replaced with something else
         /From/Size                    # voxel dimensions
         /From/Scales                  # voxel pixdims
-        /From/Mapping/Type            # "linear" - could be also be "nonlinear"
+        /From/Mapping/Type            # "linear"
         /From/Mapping/Transform       # source voxel-to-world sform
         /From/Mapping/Inverse         # optional inverse
 
@@ -121,7 +121,7 @@ def writeLinearX5(fname, xform, src, ref):
         /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)
@@ -160,24 +160,37 @@ def writeNonLinearX5(fname, field):
         /Version                      # "0.0.1"
         /Metadata                     # json string containing unstructured metadata
 
-        /Type                         # "nonlinear"
         /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/Type            # "linear"
         /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/Type              # "linear"
         /To/Mapping/Transform         # reference voxel-to-world sform
+        /To/Mapping/Type              # "linear"
         /To/Mapping/Inverse           # optional inverse
-    """
+    """  # noqa
 
     # TODO coefficient fields