From ca204fed5e767f3101c2becd4431c9c29c58a48c Mon Sep 17 00:00:00 2001 From: Paul McCarthy <pauldmccarthy@gmail.com> Date: Thu, 18 Jul 2019 11:47:02 +0100 Subject: [PATCH] RF,BK: Change nomenclature in nonlinear module. More to come. --- fsl/transform/nonlinear.py | 331 +++++++++++++++++-------------------- 1 file changed, 155 insertions(+), 176 deletions(-) diff --git a/fsl/transform/nonlinear.py b/fsl/transform/nonlinear.py index 0ea01038b..6cc516e55 100644 --- a/fsl/transform/nonlinear.py +++ b/fsl/transform/nonlinear.py @@ -8,17 +8,18 @@ FNIRT-style nonlinear transformations. -The :class:`DisplacementField` and :class:`CoefficientField` can be used to +The :class:`DeformationField` and :class:`CoefficientField` can be used to load and interact with FNIRT transformation images. The following utility functions are also available: + .. autosummary:: :nosignatures: - detectDisplacementType - convertDisplacementType - convertDisplacementSpace - coefficientFieldToDisplacementField + detectDeformationType + convertDeformationType + convertDeformationSpace + coefficientFieldToDeformationField """ @@ -38,7 +39,7 @@ log = logging.getLogger(__name__) class NonLinearTransform(fslimage.Image): """Class which represents a nonlinear transformation. This is just a base - class for the :class:`DisplacementField` and :class:`CoefficientField` + class for the :class:`DeformationField` and :class:`CoefficientField` classes. @@ -53,17 +54,6 @@ class NonLinearTransform(fslimage.Image): the corresponding location in the source image space. Therefore, these non-linear transformation effectively encode a transformation *from* the reference image *to* the source image. - - - A FNIRT nonlinear transformation often contains a *premat*, a global - affine transformation from the source space to the reference space, which - was calculated with FLIRT, and used as the starting point for the - non-linear optimisation performed by FNIRT. - - - This affine may be provided when creating a ``NonLinearTransform`` as the - ``srcToRefMat`` argument to :meth:`__init__`, and is subsequently accessed - via the :meth:`srcToRefMat` attribute. """ @@ -73,7 +63,7 @@ class NonLinearTransform(fslimage.Image): ref, srcSpace=None, refSpace=None, - srcToRefMat=None, + **kwargs): """Create a ``NonLinearTransform``. @@ -93,20 +83,11 @@ class NonLinearTransform(fslimage.Image): ``NonLinearTransform`` maps from. Defaults to ``'fsl'``. - :arg srcToRefMat: Optional initial global affine transformation from - 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 ``refSpace`` coordinates - (typically ``'fsl'`` coordinates, i.e. scaled voxels - potentially with a left-right flip). - All other arguments are passed through to :meth:`.Image.__init__`. """ - if srcSpace is None: srcSpace = 'fsl' - if refSpace is None: refSpace = 'fsl' - if srcToRefMat is not None: srcToRefMat = np.copy(srcToRefMat) + if srcSpace is None: srcSpace = 'fsl' + if refSpace is None: refSpace = 'fsl' if srcSpace not in ('fsl', 'voxel', 'world') or \ refSpace not in ('fsl', 'voxel', 'world'): @@ -119,11 +100,6 @@ class NonLinearTransform(fslimage.Image): self.__ref = fslimage.Nifti(ref.header.copy()) self.__srcSpace = srcSpace self.__refSpace = refSpace - self.__srcToRefMat = srcToRefMat - self.__refToSrcMat = None - - if srcToRefMat is not None: - self.__refToSrcMat = affine.invert(srcToRefMat) @property @@ -158,23 +134,7 @@ class NonLinearTransform(fslimage.Image): return self.__refSpace - @property - def srcToRefMat(self): - """Return the initial source-to-reference affine, or ``None`` if - there isn't one. - """ - return self.__srcToRefMat - - - @property - def refToSrcMat(self): - """Return the inverse of the initial source-to-reference affine, or - ``None`` if there isn't one. - """ - return self.__refToSrcMat - - - def transform(self, coords, from_=None, to=None, premat=True): + def transform(self, coords, from_=None, to=None): """Transform coordinates from the reference image space to the source image space. Implemented by sub-classes. @@ -186,31 +146,31 @@ class NonLinearTransform(fslimage.Image): :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() -class DisplacementField(NonLinearTransform): - """Class which represents a displacement field which, at each voxel, - contains an absolute or relative displacement between a source space and a - reference space. +class DeformationField(NonLinearTransform): + """Class which represents a deformation (a.k.a. warp) field which, at each + voxel, contains either: + + - a relative displacement from the reference space to the source space, + or + - absolute coordinates in the source space """ def __init__(self, image, src, ref=None, **kwargs): """Create a ``DisplacementField``. - :arg ref: Optional. If not provided, it is assumed that the - reference is defined in the same space as ``image``. + :arg ref: Optional. If not provided, it is assumed that the + reference is defined in the same space as ``image``. - :arg dispType: Either ``'absolute'`` or ``'relative'``, indicating - the type of this displacement field. If not provided, - will be inferred via the :func:`detectDisplacementType` - function. + :arg defType: Either ``'absolute'`` or ``'relative'``, indicating + the type of this displacement field. If not provided, + will be inferred via the :func:`detectDeformationType` + function. All other arguments are passed through to :meth:`NonLinearTransform.__init__`. @@ -219,46 +179,46 @@ class DisplacementField(NonLinearTransform): if ref is None: ref = self - dispType = kwargs.pop('dispType', None) + defType = kwargs.pop('defType', None) - if dispType not in (None, 'relative', 'absolute'): - raise ValueError('Invalid value for dispType: {}'.format(dispType)) + if defType not in (None, 'relative', 'absolute'): + raise ValueError('Invalid value for defType: {}'.format(defType)) NonLinearTransform.__init__(self, image, src, ref, **kwargs) if not self.sameSpace(self.ref): raise ValueError('Invalid reference image: {}'.format(self.ref)) - self.__dispType = dispType + self.__defType = defType @property - def displacementType(self): - """The type of this ``DisplacementField`` - ``'absolute'`` or + def deformationType(self): + """The type of this ``DeformationField`` - ``'absolute'`` or ``'relative'``. """ - if self.__dispType is None: - self.__dispType = detectDisplacementType(self) - return self.__dispType + if self.__defType is None: + self.__defType = detectDeformationType(self) + return self.__defType @property def absolute(self): - """``True`` if this ``DisplacementField`` contains absolute - displacements. + """``True`` if this ``DeformationField`` contains absolute + coordinates. """ - return self.displacementType == 'absolute' + return self.deformationType == 'absolute' @property def relative(self): - """``True`` if this ``DisplacementField`` contains relative + """``True`` if this ``DeformationField`` contains relative displacements. """ - return self.displacementType == 'relative' + return self.deformationType == 'relative' - def transform(self, coords, from_=None, to=None, premat=True): + def transform(self, coords, from_=None, to=None): """Transform the given XYZ coordinates from the reference image space to the source image space. @@ -270,9 +230,6 @@ class DisplacementField(NonLinearTransform): :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 """ @@ -310,28 +267,14 @@ class DisplacementField(NonLinearTransform): if self.absolute: disps = self.data[xs, ys, zs, :] else: disps = self.data[xs, ys, zs, :] + coords[voxmask] - # Make sure the coordinates are - # in the requested source image - # space, and apply the initial - # inv(srcToRefMat) if it there - # is one - - # And here the premat becomes - # a postmat; how confusing - if premat: postmat = self.refToSrcMat - else: postmat = None - + # Make sure the coordinates are in + # the requested source image space 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, postmat) + disps = affine.transform(disps, xform) - # Nans for input coordinates - # which were outside of the - # field + # Nans for input coordinates which + # were outside of the field outcoords = np.full(coords.shape, np.nan) outcoords[voxmask] = disps @@ -343,6 +286,17 @@ class CoefficientField(NonLinearTransform): The :meth:`displacements` method can be used to calculate relative displacements for a set of reference space voxel coordinates. + + + A FNIRT nonlinear transformation often contains a *premat*, a global + affine transformation from the source space to the reference space, which + was calculated with FLIRT, and used as the starting point for the + non-linear optimisation performed by FNIRT. + + + This affine may be provided when creating a ``CoefficientField`` as the + ``srcToRefMat`` argument to :meth:`__init__`, and is subsequently accessed + via the :meth:`srcToRefMat` attribute. """ @@ -355,6 +309,7 @@ class CoefficientField(NonLinearTransform): fieldType, knotSpacing, fieldToRefMat, + srcToRefMat=None, **kwargs): """Create a ``CoefficientField``. @@ -367,6 +322,14 @@ class CoefficientField(NonLinearTransform): image voxel coordinates into coefficient field voxel coordinates. + :arg srcToRefMat: Optional initial global affine transformation from + 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 ``refSpace`` + coordinates (typically ``'fsl'`` coordinates, i.e. + scaled voxels potentially with a left-right flip). + See the :class:`NonLinearTransform` class for details on the other arguments. """ @@ -374,6 +337,9 @@ class CoefficientField(NonLinearTransform): if fieldType not in ('cubic',): raise ValueError('Unsupported field type: {}'.format(fieldType)) + if srcToRefMat is not None: + srcToRefMat = np.copy(srcToRefMat) + NonLinearTransform.__init__(self, image, src, @@ -384,9 +350,14 @@ class CoefficientField(NonLinearTransform): self.__fieldType = fieldType self.__knotSpacing = tuple(knotSpacing) + self.__refToSrcMat = None + self.__srcToRefMat = srcToRefMat self.__fieldToRefMat = np.copy(fieldToRefMat) self.__refToFieldMat = affine.invert(self.__fieldToRefMat) + if srcToRefMat is not None: + self.__refToSrcMat = affine.invert(srcToRefMat) + @property def fieldType(self): @@ -420,11 +391,27 @@ class CoefficientField(NonLinearTransform): return np.copy(self.__refToFieldMat) + @property + def srcToRefMat(self): + """Return the initial source-to-reference affine, or ``None`` if + there isn't one. + """ + return self.__srcToRefMat + + + @property + def refToSrcMat(self): + """Return the inverse of the initial source-to-reference affine, or + ``None`` if there isn't one. + """ + return self.__refToSrcMat + + @memoize.Instanceify(memoize.memoize) - def asDisplacementField(self, dispType='relative', premat=True): - """Convert this ``CoefficientField`` to a :class:`DisplacementField`. + def asDeformationField(self, defType='relative', premat=True): + """Convert this ``CoefficientField`` to a :class:`DeformationField`. """ - return coefficientFieldToDisplacementField(self, dispType, premat) + return coefficientFieldToDeformationField(self, defType, premat) def transform(self, coords, from_=None, to=None, premat=True): @@ -446,7 +433,7 @@ class CoefficientField(NonLinearTransform): :returns: ``coords``, transformed into the source image space """ - df = self.asDisplacementField(premat=premat) + df = self.asDeformationField(premat=premat) return df.transform(coords, from_, to) @@ -531,23 +518,23 @@ class CoefficientField(NonLinearTransform): return disps -def detectDisplacementType(field): - """Attempt to automatically determine whether a displacement field is +def detectDeformationType(field): + """Attempt to automatically determine whether a deformation field is specified in absolute or relative coordinates. - :arg field: A :class:`DisplacementField` + :arg field: A :class:`DeformationField` :returns: ``'absolute'`` if it looks like ``field`` contains absolute - displacements, ``'relative'`` otherwise. + coordinates, ``'relative'`` otherwise. """ # This test is based on the assumption - # that a displacement field containing + # that a deformation field containing # absolute coordinates will have a # greater standard deviation than one # which contains relative coordinates. absdata = field[:] - reldata = convertDisplacementType(field, 'relative') + reldata = convertDeformationType(field, 'relative') stdabs = absdata.std(axis=(0, 1, 2)).sum() stdrel = reldata.std(axis=(0, 1, 2)).sum() @@ -555,20 +542,19 @@ def detectDisplacementType(field): else: return 'relative' -def convertDisplacementType(field, dispType=None): - """Convert a displacement field between storing absolute and relative - displacements. +def convertDeformationType(field, defType=None): + """Convert a deformation field between storing absolute coordinates or + relative displacements. - :arg field: A :class:`DisplacementField` instance - :arg dispType: Either ``'absolute'`` or ``'relative'``. If not provided, - the opposite type to ``field.displacementType`` is used. - :returns: A ``numpy.array`` containing the adjusted displacement - field. + :arg field: A :class:`DeformationField` instance + :arg defType: Either ``'absolute'`` or ``'relative'``. If not provided, + the opposite type to ``field.deformationType`` is used. + :returns: A ``numpy.array`` containing the adjusted deformation field. """ - if dispType is None: - if field.displacementType == 'absolute': dispType = 'relative' - else: dispType = 'absolute' + if defType is None: + if field.deformationType == 'absolute': defType = 'relative' + else: defType = 'absolute' # Regardless of the conversion direction, # we need the coordinates of every voxel @@ -588,36 +574,29 @@ def convertDisplacementType(field, dispType=None): # (what is assumed to be) the relative shift. # Or, to convert from absolute to relative, # we subtract the reference image voxels. - if dispType == 'absolute': return field.data + coords - elif dispType == 'relative': return field.data - coords + if defType == 'absolute': return field.data + coords + elif defType == 'relative': return field.data - coords -def convertDisplacementSpace(field, from_, to): - """Adjust the source and/or reference spaces of the given displacement +def convertDeformationSpace(field, from_, to): + """Adjust the source and/or reference spaces of the given deformation field. See the :meth:`.Nifti.getAffine` method for the valid values for the ``from_`` and ``to`` arguments. - :arg field: A :class:`DisplacementField` instance + :arg field: A :class:`DeformationField` instance :arg from_: New reference image coordinate system :arg to: New source image coordinate system - :returns: A new :class:`DisplacementField` which transforms between + :returns: A new :class:`DeformationField` which transforms between the reference ``from_`` coordinate system and the source ``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. fieldcoords = field.data - if field.relative: srccoords = convertDisplacementType(field) + if field.relative: srccoords = convertDeformationType(field) else: srccoords = fieldcoords srccoords = srccoords.reshape((-1, 3)) @@ -631,13 +610,13 @@ def convertDisplacementSpace(field, from_, to): srccoords = affine.transform(srccoords, srcmat) # If we have been asked to return - # an absolute displacement, the + # absolute coordinates, the # reference "from_" coordinate # system is irrelevant - we're done. if field.absolute: fieldcoords = srccoords - # Otherwise our displacement field + # Otherwise our deformation field # will contain relative displacements # between the reference image "from_" # coordinate system and the source @@ -659,30 +638,30 @@ def convertDisplacementSpace(field, from_, to): fieldcoords = srccoords - refcoords - return DisplacementField( + return DeformationField( fieldcoords.reshape(field.shape), header=field.header, src=field.src, ref=field.ref, srcSpace=to, refSpace=from_, - dispType=field.displacementType) + defType=field.deformationType) -def coefficientFieldToDisplacementField(field, - dispType='relative', - premat=True): - """Convert a :class:`CoefficientField` into a :class:`DisplacementField`. +def coefficientFieldToDeformationField(field, + defType='relative', + premat=True): + """Convert a :class:`CoefficientField` into a :class:`DeformationField`. - :arg field: :class:`CoefficientField` to convert + :arg field: :class:`CoefficientField` to convert - :arg dispType: The type of displacement field - either ``'relative'`` (the - default) or ``'absolute'``. + :arg defType: The type of deformation field - either ``'relative'`` (the + default) or ``'absolute'``. - :arg premat: If ``True`` (the default), the :meth:`srcToRefMat` is - encoded into the displacements. + :arg premat: If ``True`` (the default), the :meth:`srcToRefMat` is + encoded into the deformation field. - :return: :class:`DisplacementField` calculated from ``field``. + :return: :class:`DeformationField` calculated from ``field``. """ # Generate coordinates for every @@ -715,22 +694,22 @@ def coefficientFieldToDisplacementField(field, # from ref space to aligned-src # space. disps = field.displacements(xyz).reshape((ix, iy, iz, 3)) - rdfield = DisplacementField(disps, - src=field.src, - ref=field.ref, - srcSpace=field.srcSpace, - refSpace=field.refSpace, - header=field.ref.header, - dispType='relative') - - if (dispType == 'relative') and (not premat): + rdfield = DeformationField(disps, + src=field.src, + ref=field.ref, + srcSpace=field.srcSpace, + refSpace=field.refSpace, + header=field.ref.header, + defType='relative') + + if (defType == 'relative') and (not premat): return rdfield # Convert to absolute - the - # displacements will now be + # deformations will now be # absolute coordinates in # aligned-src space - disps = convertDisplacementType(rdfield) + disps = convertDeformationType(rdfield) # Apply the premat if requested - # this will transform the coordinates @@ -756,23 +735,23 @@ def coefficientFieldToDisplacementField(field, # # disps = affine.transform(disps, refToSrc) - adfield = DisplacementField(disps, - src=field.src, - ref=field.ref, - srcSpace=field.srcSpace, - refSpace=field.refSpace, - header=field.ref.header, - dispType='absolute') + adfield = DeformationField(disps, + src=field.src, + ref=field.ref, + srcSpace=field.srcSpace, + refSpace=field.refSpace, + header=field.ref.header, + defType='absolute') # Not either return absolute displacements, # or convert back to relative displacements - if dispType == 'absolute': + if defType == 'absolute': return adfield else: - return DisplacementField(convertDisplacementType(adfield), - src=field.src, - ref=field.ref, - srcSpace=field.srcSpace, - refSpace=field.refSpace, - header=field.ref.header, - dispType='relative') + return DeformationField(convertDeformationType(adfield), + src=field.src, + ref=field.ref, + srcSpace=field.srcSpace, + refSpace=field.refSpace, + header=field.ref.header, + defType='relative') -- GitLab