RF: toFnirt will not convert coefficient fields to displacement fields if it

can get away with it. Can't be bothered figuring out how to adjust srcToRefMat
in convertDisplacementSpace, so am just not supporting it

fsl/transform/fnirt.py

@@ -18,6 +18,9 @@ transformation matrices. The following functions are available:
 
 import logging
 
+import numpy   as np
+import nibabel as nib
+
 import fsl.data.constants as constants
 
 
@@ -154,25 +157,106 @@ def readFnirt(fname, src, ref, dispType=None):
 
 def toFnirt(field):
     """Convert a :class:`.NonLinearTransform` to a FNIRT-compatible
-    :class:`.DisplacementField`.
+    :class:`.DisplacementField` or:class:`.CoefficientField`.
 
     :arg field: :class:`.NonLinearTransform` to convert
-    :return:    A FNIRT-compatible :class:`.DisplacementField`.
+    :return:    A FNIRT-compatible :class:`.DisplacementField` or
+                :class:`.CoefficientField`.
     """
 
     from . import nonlinear
 
-    # 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)
+    # 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 displacement field.
+    #
+    # 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.
+    else:
 
-    field = nonlinear.convertDisplacementSpace(field, from_='fsl', to='fsl')
-    field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
+        if isinstance(field, nonlinear.CoefficientField):
+            field = nonlinear.coefficientFieldToDisplacementField(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.DisplacementField(
+                field.data,
+                src=field.src,
+                ref=field.ref,
+                fieldType=field.fieldType,
+                dispType=field.displacementType)
+
+        # Otherwise we have to adjust the
+        # displacements so they transform
+        # between fsl coordinates.
+        field = nonlinear.convertDisplacementSpace(
+            field, from_='fsl', to='fsl')
+
+        field.header['intent_code'] = constants.FSL_FNIRT_DISPLACEMENT_FIELD
 
     return field
 

fsl/transform/nonlinear.py

@@ -97,7 +97,7 @@ class NonLinearTransform(fslimage.Image):
                           the source image to the reference image. This is
                           assumed to be a FLIRT-style matrix, i.e. it
                           transforms from source image ``srcSpace`` coordinates
-                          into reference image ``srcSpace`` coordinates
+                          into reference image ``refSpace`` coordinates
                           (typically ``'fsl'`` coordinates, i.e. scaled voxels
                           potentially with a left-right flip).
 
@@ -174,15 +174,21 @@ class NonLinearTransform(fslimage.Image):
         return self.__refToSrcMat
 
 
-    def transform(self, coords, from_=None, to=None):
+    def transform(self, coords, from_=None, to=None, premat=True):
         """Transform coordinates from the reference image space to the source
         image space. Implemented by sub-classes.
 
         :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
                      ``(n, 3)`` containing ``n`` sets of coordinates in the
                      reference space.
+
         :arg from_:  Reference image space that ``coords`` are defined in
+
         :arg to:     Source image space to transform ``coords`` into
+
+        :arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
+                     to the coordinates after they have been calculated.
+
         :returns:    ``coords``, transformed into the source image space
         """
         raise NotImplementedError()
@@ -252,15 +258,21 @@ class DisplacementField(NonLinearTransform):
         return self.displacementType == 'relative'
 
 
-    def transform(self, coords, from_=None, to=None):
+    def transform(self, coords, from_=None, to=None, premat=True):
         """Transform the given XYZ coordinates from the reference image space
         to the source image space.
 
         :arg coords: A sequence of XYZ coordinates, or ``numpy`` array of shape
                      ``(n, 3)`` containing ``n`` sets of coordinates in the
                      reference space.
+
         :arg from_:  Reference image space that ``coords`` are defined in
+
         :arg to:     Source image space to transform ``coords`` into
+
+        :arg premat: If ``True``, the inverse :meth:`srcToRefMat` is applied
+                     to the coordinates after they have been calculated.
+
         :returns:    ``coords``, transformed into the source image space
         """
 
@@ -303,14 +315,19 @@ class DisplacementField(NonLinearTransform):
         # space, and apply the initial
         # inv(srcToRefMat) if it there
         # is one
-        postmat = self.refToSrcMat
+
+        # And here the premat becomes
+        # a postmat; how confusing
+        if premat: postmat = self.refToSrcMat
+        else:      postmat = None
+
         if to != self.srcSpace:
             xform = self.src.getAffine(self.srcSpace, to)
             if postmat is not None: postmat = affine.concat(xform, postmat)
             else:                   postmat = xform
 
         if postmat is not None:
-            disps = affine.transform(disps, xform)
+            disps = affine.transform(disps, postmat)
 
         # Nans for input coordinates
         # which were outside of the
@@ -322,8 +339,7 @@ class DisplacementField(NonLinearTransform):
 
 
 class CoefficientField(NonLinearTransform):
-    """Class which represents a cubic B-spline coefficient field generated by
-    FNIRT.
+    """Class which represents a B-spline coefficient field generated by FNIRT.
 
     The :meth:`displacements` method can be used to calculate relative
     displacements for a set of reference space voxel coordinates.
@@ -590,6 +606,13 @@ def convertDisplacementSpace(field, from_, to):
                 coordinate system.
     """
 
+    # We can't adjust the displacements
+    # for fields with an initial
+    # source-to-ref affine
+    if field.srcToRefMat is not None:
+        raise ValueError('Cannot adjust displacements of '
+                         'fields with an initial affine')
+
     # Get the field in absolute coordinates
     # if necessary - these are our source
     # coordinates in the original "to" space.