diff --git a/fsl/transform/fnirt.py b/fsl/transform/fnirt.py
index 79c0edbb9a7acb58b376f09eff91bde289dab1b1..30db64fd34e7d82dbc915859025388ff97822628 100644
--- a/fsl/transform/fnirt.py
+++ b/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
diff --git a/fsl/transform/nonlinear.py b/fsl/transform/nonlinear.py
index a32225de23e18e031447a27e79ded8c552cd19d4..0ea01038beed4e17589a1f35e0fdfd7a76beaac9 100644
--- a/fsl/transform/nonlinear.py
+++ b/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.