From e800b3d158875b77dc13404397f04d6891a34bbf Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Sat, 13 Jul 2019 17:00:12 +0100
Subject: [PATCH] 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 | 108 ++++++++++++++++++++++++++++++++-----
fsl/transform/nonlinear.py | 37 ++++++++++---
2 files changed, 126 insertions(+), 19 deletions(-)
diff --git a/fsl/transform/fnirt.py b/fsl/transform/fnirt.py
index 79c0edbb9..30db64fd3 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 a32225de2..0ea01038b 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.
--
GitLab