fsl/transform/affine.py

@@ -13,6 +13,7 @@ transformations. The following functions are available:
    transform
    scaleOffsetXform
    invert
+   flip
    concat
    compose
    decompose
@@ -20,6 +21,7 @@ transformations. The following functions are available:
    rotMatToAxisAngles
    axisAnglesToRotMat
    axisBounds
+   mergeBounds
    rmsdev
    rescale
 
@@ -60,7 +62,7 @@ def veclength(vec):
 
     Multiple vectors may be passed in, with a shape of ``(n, 3)``.
     """
-    vec = np.array(vec, copy=False).reshape(-1, 3)
+    vec = np.asarray(vec).reshape(-1, 3)
     return np.sqrt(np.einsum('ij,ij->i', vec, vec))
 
 
@@ -69,7 +71,7 @@ def normalise(vec):
 
     Multiple vectors may be passed in, with a shape of ``(n, 3)``.
     """
-    vec = np.array(vec, copy=False).reshape(-1, 3)
+    vec = np.asarray(vec).reshape(-1, 3)
     n   = (vec.T / veclength(vec)).T
 
     if n.size == 3:
@@ -78,7 +80,41 @@ def normalise(vec):
     return n
 
 
-def scaleOffsetXform(scales, offsets):
+def flip(shape, xform, *axes):
+    """Applies a flip/inversion to an affine transform along specified axes.
+
+    :arg shape: The ``(x, y, z)`` shape of the data.
+    :arg xform: Transformation matrix which transforms voxel coordinates
+                to a world coordinate system.
+    :arg axes:  Indices of axes to flip
+    """
+
+    if len(axes) == 0:
+        return xform
+
+    axes = sorted(set(axes))
+    if any(a not in (0, 1, 2) for a in axes):
+        raise ValueError(f'Invalid axis: {axes}')
+
+    los, his = axisBounds(shape, xform, axes, boundary=None)
+    scales   = [1, 1, 1]
+    pre      = [0, 0, 0]
+    post     = [0, 0, 0]
+
+    for ax, lo, hi in zip(axes, los, his):
+        mid        = lo + (hi - lo) / 2
+        scales[ax] = -1
+        pre[   ax] = -mid
+        post[  ax] =  mid
+
+    pre    = scaleOffsetXform(None, pre)
+    post   = scaleOffsetXform(None, post)
+    scales = scaleOffsetXform(scales)
+
+    return concat(post, scales, pre, xform)
+
+
+def scaleOffsetXform(scales=None, offsets=None):
     """Creates and returns an affine transformation matrix which encodes
     the specified scale(s) and offset(s).
 
@@ -94,6 +130,9 @@ def scaleOffsetXform(scales, offsets):
     :returns:     A ``numpy.float32`` array of size :math:`4 \\times 4`.
     """
 
+    if scales  is None: scales  = [1, 1, 1]
+    if offsets is None: offsets = [0, 0, 0]
+
     oktypes = (abc.Sequence, np.ndarray)
 
     if not isinstance(scales,  oktypes): scales  = [scales]
@@ -120,19 +159,28 @@ def scaleOffsetXform(scales, offsets):
     return xform
 
 
-def compose(scales, offsets, rotations, origin=None):
+def compose(scales, offsets, rotations, origin=None, shears=None,
+            scaleAtOrigin=False):
     """Compose a transformation matrix out of the given scales, offsets
     and axis rotations.
 
-    :arg scales:    Sequence of three scale values.
+    :arg scales:        Sequence of three scale values.
+
+    :arg offsets:       Sequence of three offset values.
 
-    :arg offsets:   Sequence of three offset values.
+    :arg rotations:     Sequence of three rotation values, in radians, or
+                        a rotation matrix of shape ``(3, 3)``.
 
-    :arg rotations: Sequence of three rotation values, in radians, or
-                    a rotation matrix of shape ``(3, 3)``.
+    :arg origin:        Origin of rotation.  If not provided, the rotation
+                        origin is ``(0, 0, 0)``.
 
-    :arg origin:    Origin of rotation - must be scaled by the ``scales``.
-                    If not provided, the rotation origin is ``(0, 0, 0)``.
+    :arg shears:        Sequence of three shear values
+
+    :arg scaleAtOrigin: If ``True``, the scaling parameters are applied with
+                        respect to the ``origin``, i.e. so that the location
+                        of ``origin`` is unchanged after scaling.  If
+                        ``False`` (the default), the scaling parameters are
+                        applied with respect to location ``(0, 0, 0)``.
     """
 
     preRotate  = np.eye(4)
@@ -154,6 +202,7 @@ def compose(scales, offsets, rotations, origin=None):
     scale  = np.eye(4, dtype=np.float64)
     offset = np.eye(4, dtype=np.float64)
     rotate = np.eye(4, dtype=np.float64)
+    shear  = np.eye(4, dtype=np.float64)
 
     scale[  0,  0] = scales[ 0]
     scale[  1,  1] = scales[ 1]
@@ -164,10 +213,18 @@ def compose(scales, offsets, rotations, origin=None):
 
     rotate[:3, :3] = rotations
 
-    return concat(offset, postRotate, rotate, preRotate, scale)
+    if shears is not None:
+        shear[0, 1] = shears[0]
+        shear[0, 2] = shears[1]
+        shear[1, 2] = shears[2]
+
+    if scaleAtOrigin:
+        return concat(offset, postRotate, rotate, scale, preRotate, shear)
+    else:
+        return concat(offset, postRotate, rotate, preRotate, scale, shear)
 
 
-def decompose(xform, angles=True):
+def decompose(xform, angles=True, shears=False):
     """Decomposes the given transformation matrix into separate offsets,
     scales, and rotations, according to the algorithm described in:
 
@@ -175,8 +232,7 @@ def decompose(xform, angles=True):
     320-323 in *Graphics Gems II*, James Arvo (editor), Academic Press, 1991,
     ISBN: 0120644819.
 
-    It is assumed that the given transform has no perspective components. Any
-    shears in the affine are discarded.
+    It is assumed that the given transform has no perspective components.
 
     :arg xform:  A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
 
@@ -184,6 +240,8 @@ def decompose(xform, angles=True):
                  as axis-angles, in radians. Otherwise, the rotation matrix
                  is returned.
 
+    :arg shears: Defaults to ``False``. If ``True``, shears are returned.
+
     :returns: The following:
 
                - A sequence of three scales
@@ -191,6 +249,7 @@ def decompose(xform, angles=True):
                  was a ``(3, 3)`` matrix)
                - A sequence of three rotations, in radians. Or, if
                  ``angles is False``, a rotation matrix.
+               - If ``shears is True``, a sequence of three shears.
     """
 
     # The inline comments in the code below are taken verbatim from
@@ -199,7 +258,7 @@ def decompose(xform, angles=True):
     # The next step is to extract the translations. This is trivial;
     # we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
     # point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
-    xform = xform.T
+    xform = np.array(xform).T
 
     if xform.shape == (4, 4):
         translations = xform[ 3, :3]
@@ -214,6 +273,7 @@ def decompose(xform, angles=True):
     # The process of finding the scaling factors and shear parameters
     # is interleaved. First, find s_x = |M'_1|.
     sx = np.sqrt(np.dot(M1, M1))
+    M1 = M1 / sx
 
     # Then, compute an initial value for the xy shear factor,
     # s_xy = M'_1 * M'_2. (this is too large by the y scaling factor).
@@ -230,7 +290,7 @@ def decompose(xform, angles=True):
     # The second row is normalized, and s_xy is divided by s_y to
     # get its final value.
     M2  = M2  / sy
-    sxy = sxy / sy
+    sxy = sxy / sx
 
     # The xz and yz shear factors are computed as in the preceding,
     sxz = np.dot(M1, M3)
@@ -245,8 +305,8 @@ def decompose(xform, angles=True):
     # the third row is normalized, and the xz and yz shear factors are
     # rescaled.
     M3  = M3  / sz
-    sxz = sxz / sz
-    syz = sxz / sz
+    sxz = sxz / sx
+    syz = syz / sy
 
     # The resulting matrix now is a pure rotation matrix, except that it
     # might still include a scale factor of -1. If the determinant of the
@@ -266,7 +326,13 @@ def decompose(xform, angles=True):
     if angles: rotations = rotMatToAxisAngles(R.T)
     else:      rotations = R.T
 
-    return np.array([sx, sy, sz]), translations, rotations
+    retval = [np.array([sx, sy, sz]), translations, rotations]
+
+    if shears:
+        retval.append(np.array((sxy, sxz, syz)))
+
+    return tuple(retval)
+
 
 
 def rotMatToAffine(rotmat, origin=None):
@@ -451,11 +517,45 @@ def axisBounds(shape,
     else:      return (lo,    hi)
 
 
+def mergeBounds(*bounds):
+    """Merge multiple ``(lo, hi)`` bounds together. Returns the union of
+    all bounds, i.e. a new set of bounds which encompasses all given bounds.
+
+    Each set of bounds can be provided in one of the following formats:
+
+      - ``((lo, lo, lo), (hi, hi, hi))``
+      - ``((lo, hi), (lo, hi), (lo, hi))``
+
+    The return value will be a tuple of tuples, otherwise having the same
+    format as the inputs.
+    """
+
+    nbounds = len(bounds)
+    bounds  = np.array(bounds)
+
+    if bounds.shape not in ((nbounds, 3, 2), (nbounds, 2, 3)):
+        raise ValueError('Unsupported bounds format')
+
+    if bounds.shape == (nbounds, 2, 3):
+        return ((bounds[:, 0, 0].min(),
+                 bounds[:, 0, 1].min(),
+                 bounds[:, 0, 2].min()),
+                (bounds[:, 1, 0].max(),
+                 bounds[:, 1, 1].max(),
+                 bounds[:, 1, 2].max()))
+    else:
+        return ((bounds[:, 0, 0].min(),
+                 bounds[:, 0, 1].max()),
+                (bounds[:, 1, 0].min(),
+                 bounds[:, 1, 1].max()),
+                (bounds[:, 2, 0].min(),
+                 bounds[:, 2, 1].max()))
+
+
 def transform(p, xform, axes=None, vector=False):
     """Transforms the given set of points ``p`` according to the given affine
     transformation ``xform``.
 
-
     :arg p:      A sequence or array of points of shape :math:`N \\times  3`.
 
     :arg xform:  A ``(4, 4)`` affine transformation matrix with which to
@@ -610,8 +710,8 @@ def rescale(oldShape, newShape, origin=None):
     if origin is None:
         origin = 'centre'
 
-    oldShape = np.array(oldShape, dtype=np.float)
-    newShape = np.array(newShape, dtype=np.float)
+    oldShape = np.array(oldShape, dtype=float)
+    newShape = np.array(newShape, dtype=float)
     ndim     = len(oldShape)
 
     if len(oldShape) != len(newShape):

fsl/transform/flirt.py

@@ -63,7 +63,7 @@ def fromFlirt(xform, src, ref, from_='voxel', to='world'):
     Valid values for the ``from_`` and ``to`` arguments are:
 
      - ``voxel``: The voxel coordinate system
-     - ``fsl``: The FSL coordiante system (voxels scaled by pixdims, with the
+     - ``fsl``: The FSL coordinate system (voxels scaled by pixdims, with the
        X axis inverted if the image sform/qform has a positive determinant)
      - ``world``  The world coordinate system
 
@@ -138,15 +138,15 @@ def sformToFlirtMatrix(srcImage, refImage, srcXform=None):
                    :attr:`.Nifti.voxToWorldMat`
     """
 
-    srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat
-    srcVoxToWorldMat     = srcImage.voxToWorldMat
-    refWorldToVoxMat     = refImage.worldToVoxMat
-    refVoxToScaledVoxMat = refImage.voxToScaledVoxMat
+    srcFSLToVoxMat   = srcImage.getAffine('fsl',   'voxel')
+    srcVoxToWorldMat = srcImage.getAffine('voxel', 'world')
+    refWorldToVoxMat = refImage.getAffine('world', 'voxel')
+    refVoxToFSLMat   = refImage.getAffine('voxel', 'fsl')
 
     if srcXform is not None:
         srcVoxToWorldMat = srcXform
 
-    return concat(refVoxToScaledVoxMat,
+    return concat(refVoxToFSLMat,
                   refWorldToVoxMat,
                   srcVoxToWorldMat,
-                  srcScaledVoxToVoxMat)
+                  srcFSLToVoxMat)

fsl/transform/fnirt.py

@@ -82,6 +82,12 @@ the reference image, and may contain:
    the source image coordinates which correspond to those reference image
    coordinates.
 
+.. note:: FNIRT deformation field files give no indication as to whether they
+          contain relative displacements or absolute coordinates, so heuristics
+          must be used to infer what is stored in a a particular file. The
+          :func:`.nonlinear.detectDeformationType` function can be used to
+          determine whether a file contains relative displacements or absolute
+          coordinates.
 
 If an initial linear registration was used as the starting point for FNIRT,
 this is encoded into the displacements/coordinates themselves, so they can be
@@ -219,6 +225,27 @@ def _readFnirtCoefficientField(fname, img, src, ref):
     # image voxel coordinates.
     fieldToRefMat = affine.scaleOffsetXform(knotSpacing, 0)
 
+    # But if the provided reference has
+    # different resolution to the
+    # reference that was originally
+    # used to generate the warp field,
+    # we need to adjust the field
+    # accordingly. We assume that the
+    # references are aligned in the FSL
+    # coordinate system, so simply apply
+    # a scaling factor calculated by
+    # dividing the original reference
+    # pixdims by the provided reference
+    # pixdims.
+    refPixdims = np.array([img.header['intent_p1'],
+                           img.header['intent_p2'],
+                           img.header['intent_p3']])
+
+    if not np.all(np.isclose(ref.pixdim[:3], refPixdims)):
+        fieldToRefMat = affine.concat(
+            affine.scaleOffsetXform(refPixdims / ref.pixdim[:3], 0),
+            fieldToRefMat)
+
     return nonlinear.CoefficientField(fname,
                                       src,
                                       ref,
@@ -230,7 +257,7 @@ def _readFnirtCoefficientField(fname, img, src, ref):
                                       fieldToRefMat=fieldToRefMat)
 
 
-def readFnirt(fname, src, ref, defType=None):
+def readFnirt(fname, src, ref, defType=None, intent=None):
     """Reads a non-linear FNIRT transformation image, returning
     a :class:`.DeformationField` or :class:`.CoefficientField` depending
     on the file type.
@@ -241,28 +268,38 @@ def readFnirt(fname, src, ref, defType=None):
     :arg defType:  Deformation type - either ``'absolute'`` or ``'relative'``.
                    Only used if the file is a deformation field. If not
                    provided, is automatically inferred from the data.
+    :arg intent:   NIFTI intent code of ``fname``. e.g.
+                   :attr:`.constants.FSL_FNIRT_DISPLACEMENT_FIELD`. If not
+                   provided, the intent is read from the image header.
     """
 
-    # Figure out whether the file
-    # is a deformation field or
-    # a coefficient field
-
-    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 defType not in (None, 'absolute', 'relative'):
+        raise ValueError('defType must be None, "absolute" or "relative" '
+                         '(passed in as {})'.format(defType))
+
+    # Figure out whether the file is a
+    # deformation field or a coefficient
+    # field by checking the intent code.
+    # If the intent is provided, assume
+    # that the caller knows the type of
+    # the field.
+    img    = fslimage.Image(fname)
+    intent = intent or img.intent
+    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:
+    if intent in disps:
         return _readFnirtDeformationField(fname, img, src, ref, defType)
-    elif img.intent in coefs:
+    elif intent in coefs:
         return _readFnirtCoefficientField(fname, img, src, ref)
     else:
-        raise ValueError('Cannot determine type of nonlinear '
-                         'file {}'.format(fname))
+        raise ValueError('Cannot determine type of nonlinear warp field '
+                         '{} (intent code: {})'.format(fname, intent))
 
 
 def toFnirt(field):

fsl/transform/nonlinear.py

@@ -24,16 +24,17 @@ utility functions are also available:
 """
 
 
-import                                logging
-import itertools                   as it
+import                             logging
+import itertools                as it
 
-import numpy                       as np
-import scipy.ndimage.interpolation as ndinterp
+import numpy                    as np
+import scipy.ndimage            as ndimage
 
-import fsl.utils.memoize           as memoize
-import fsl.data.image              as fslimage
-import fsl.utils.image.resample    as resample
-from . import                         affine
+import fsl.utils.memoize        as memoize
+import fsl.data.image           as fslimage
+import fsl.data.constants       as fslconstants
+import fsl.utils.image.resample as resample
+from . import                      affine
 
 
 log = logging.getLogger(__name__)
@@ -200,10 +201,29 @@ class DeformationField(NonLinearTransform):
         defType = kwargs.pop('defType', None)
 
         if defType not in (None, 'relative', 'absolute'):
-            raise ValueError('Invalid value for defType: {}'.format(defType))
+            raise ValueError(f'Invalid value for defType: {defType}')
 
         NonLinearTransform.__init__(self, image, src, ref, **kwargs)
 
+        # Displacement fields generated by the FSL
+        # convertwarp tool prior to FSL 6.0.5 would
+        # often not have their sform/qform affines
+        # set correctly. Here we check the sform/
+        # qform and, if necessary, construct a new
+        # affine from the reference image.
+        if self.getXFormCode() == fslconstants.NIFTI_XFORM_UNKNOWN:
+            xform = self.ref.voxToWorldMat
+
+            # It is possible for a displacement field
+            # to have different resolution to its
+            # reference space - they just have to be
+            # aligned in the world coordinate system.
+            if self.shape[:3] != self.ref.shape[:3]:
+                rescale = affine.rescale(self.shape[:3], self.ref.shape[:3])
+                xform   = affine.concat(xform, rescale)
+
+            self.voxToWorldMat = xform
+
         self.__defType = defType
 
 
@@ -251,7 +271,9 @@ class DeformationField(NonLinearTransform):
         if from_ is None: from_ = self.refSpace
         if to    is None: to    = self.srcSpace
 
-        coords = np.asanyarray(coords)
+        coords   = np.asanyarray(coords)
+        outshape = coords.shape
+        coords   = coords.reshape((-1, 3))
 
         # We may need to pre-transform the
         # coordinates so they are in the
@@ -278,7 +300,7 @@ class DeformationField(NonLinearTransform):
         # Mask out the coordinates
         # that are out of bounds of
         # the deformation field
-        voxels  = np.round(voxels).astype(np.int)
+        voxels  = np.round(voxels).astype(np.int32)
         voxmask = (voxels >= [0, 0, 0]) & (voxels < self.shape[:3])
         voxmask = voxmask.all(axis=1)
         voxels  = voxels[voxmask]
@@ -299,7 +321,7 @@ class DeformationField(NonLinearTransform):
         outcoords          = np.full(coords.shape, np.nan)
         outcoords[voxmask] = disps
 
-        return outcoords
+        return outcoords.reshape(outshape)
 
 
 class CoefficientField(NonLinearTransform):
@@ -494,13 +516,11 @@ class CoefficientField(NonLinearTransform):
 
         fdata      = self.data
         nx, ny, nz = self.shape[:3]
-        ix, iy, iz = self.ref.shape[:3]
 
         # Convert the given voxel coordinates
         # into the corresponding coefficient
         # field voxel coordinates
-        x, y, z    = coords.T
-        i, j, k    = affine.transform(coords, self.refToFieldMat).T
+        i, j, k = affine.transform(coords, self.refToFieldMat).T
 
         # i, j, k: coefficient field indices
         # u, v, w: position of the ref voxel
@@ -508,9 +528,9 @@ class CoefficientField(NonLinearTransform):
         u = np.remainder(i, 1)
         v = np.remainder(j, 1)
         w = np.remainder(k, 1)
-        i = np.floor(i).astype(np.int)
-        j = np.floor(j).astype(np.int)
-        k = np.floor(k).astype(np.int)
+        i = np.floor(i).astype(np.int32)
+        j = np.floor(j).astype(np.int32)
+        k = np.floor(k).astype(np.int32)
 
         disps = np.zeros(coords.shape)
 
@@ -580,6 +600,10 @@ def convertDeformationType(field, defType=None):
         if field.deformationType == 'absolute': defType = 'relative'
         else:                                   defType = 'absolute'
 
+    if defType not in ('absolute', 'relative'):
+        raise ValueError('defType must be "absolute" or "relative" '
+                         '("{}" passed)'.format(defType))
+
     # Regardless of the conversion direction,
     # we need the coordinates of every voxel
     # in the reference coordinate system.
@@ -599,8 +623,8 @@ def convertDeformationType(field, defType=None):
     # assumed to be) the relative shift. Or,
     # to convert from absolute to relative,
     # we subtract the reference image voxels.
-    if   defType == 'absolute': return field.data + coords
-    elif defType == 'relative': return field.data - coords
+    if defType == 'absolute': return field.data + coords
+    else:                     return field.data - coords
 
 
 def convertDeformationSpace(field, from_, to):
@@ -679,7 +703,13 @@ def convertDeformationSpace(field, from_, to):
         defType=field.deformationType)
 
 
-def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
+def applyDeformation(image,
+                     field,
+                     ref=None,
+                     order=1,
+                     mode=None,
+                     cval=None,
+                     premat=None):
     """Applies a :class:`DeformationField` to an :class:`.Image`.
 
     The image is transformed into the space of the field's reference image
@@ -691,25 +721,30 @@ def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
     the input image. It is therefore assumed that an alternate ``ref`` is
     aligned in world coordinates with the field's actual reference image.
 
-    :arg image: :class:`.Image` to be transformed
+    :arg image:  :class:`.Image` to be transformed
 
-    :arg field: :class:`DeformationField` to use
+    :arg field:  :class:`DeformationField` to use
 
-    :arg ref:   Alternate reference image - if not provided, ``field.ref``
-                is used
+    :arg ref:    Alternate reference image - if not provided, ``field.ref``
+                 is used
 
-    :arg order: Spline interpolation order, passed through to the
-                ``scipy.ndimage.affine_transform`` function - ``0``
-                corresponds to nearest neighbour interpolation, ``1``
-                (the default) to linear interpolation, and ``3`` to
-                cubic interpolation.
+    :arg order:  Spline interpolation order, passed through to the
+                 ``scipy.ndimage.affine_transform`` function - ``0``
+                 corresponds to nearest neighbour interpolation, ``1``
+                 (the default) to linear interpolation, and ``3`` to
+                 cubic interpolation.
 
-    :arg mode:  How to handle regions which are outside of the image FOV.
-                Defaults to `''nearest'``.
+    :arg mode:   How to handle regions which are outside of the image FOV.
+                 Defaults to `''nearest'``.
 
-    :arg cval:  Constant value to use when ``mode='constant'``.
+    :arg cval:   Constant value to use when ``mode='constant'``.
 
-    :return:    ``numpy.array`` containing the transformed image data.
+    :arg premat: Optional affine transform which can be used if ``image`` is
+                 not in the same space as ``field.src``. Assumed to transform
+                 from ``image`` **voxel** coordinates into ``field.src``
+                 **voxel** coordinates.
+
+    :return:     ``numpy.array`` containing the transformed image data.
     """
 
     if order is None: order = 1
@@ -763,20 +798,23 @@ def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
     # We assume world-world alignment
     # between the original source
     # and the image to be resampled
-    if not image.sameSpace(src):
-        post  = affine.concat(image.getAffine('world', 'voxel'),
-                              src  .getAffine('voxel', 'world'))
+    if (premat is not None) or (not image.sameSpace(src)):
+        if premat is None:
+            premat = affine.concat(image.getAffine('world', 'voxel'),
+                                   src  .getAffine('voxel', 'world'))
+        else:
+            premat = affine.invert(premat)
         shape = field.shape
         field = field.reshape((-1, 3))
-        field = affine.transform(field, post)
+        field = affine.transform(field, premat)
         field = field.reshape(shape)
 
     field = field.transpose((3, 0, 1, 2))
-    return ndinterp.map_coordinates(image.data,
-                                    field,
-                                    order=order,
-                                    mode=mode,
-                                    cval=cval)
+    return ndimage.map_coordinates(image.data,
+                                   field,
+                                   order=order,
+                                   mode=mode,
+                                   cval=cval)
 
 
 def coefficientFieldToDeformationField(field, defType='relative', premat=True):

fsl/utils/assertions.py

@@ -33,25 +33,32 @@ import fsl.utils.ensure         as ensure
 import fsl.data.melodicanalysis as fslma
 
 
-_DISABLE_ASSERTIONS = False
-"""
-"""
+_DISABLE_ASSERTIONS = 0
+"""Semaphore used by the :func:`disabled` context manager. """
 
 
 @contextlib.contextmanager
-def disabled():
+def disabled(disable=True):
     """Context manager which allows assertion checks to be temporarily
     disabled.
+
+    If calls to this function are nested, only one of the calls need to be made
+    with ``disable=True`` for assertions to be disabled; any other calls which
+    are part of the call stack which set ``disable=False`` will have no effect.
+
+    :arg disable: Set to ``True`` (the default) to disable assertions,
+                  or ``False`` to enable them.
     """
     global _DISABLE_ASSERTIONS
 
-    oldval              = _DISABLE_ASSERTIONS
-    _DISABLE_ASSERTIONS = True
+    if disable:
+        _DISABLE_ASSERTIONS += 1
 
     try:
         yield
     finally:
-        _DISABLE_ASSERTIONS = oldval
+        if disable:
+            _DISABLE_ASSERTIONS -= 1
 
 
 def _canDisable(func):
@@ -59,7 +66,7 @@ def _canDisable(func):
     via the :func:`disabled` context manager.
     """
     def wrapper(*args, **kwargs):
-        if not _DISABLE_ASSERTIONS:
+        if _DISABLE_ASSERTIONS == 0:
             return func(*args, **kwargs)
     return wrapper
 

fsl/utils/bids.py

@@ -9,6 +9,7 @@
 .. autosummary::
    :nosignatures:
 
+   BIDSFile
    isBIDSDir
    inBIDSDir
    isBIDSFile
@@ -17,6 +18,8 @@
 All of the other functions in this module should not be considered part of the
 public API.
 
+See https://bids-specification.readthedocs.io/en/stable/ for more information
+about BIDS.
 
 .. note::  The `pybids <https://bids-standard.github.io/pybids/>`_ library is
            a more suitable choice if you are after a more robust and featured
@@ -35,8 +38,8 @@ import fsl.utils.path    as fslpath
 
 
 class BIDSFile(object):
-    """The ``BIDSFile`` class parses and stores the entities and suffix contained
-    in a BIDS file. See the :func:`parseFilename` function.
+    """The ``BIDSFile`` class parses and stores the entities and suffix
+    contained in a BIDS file. See the :func:`parseFilename` function.
 
     The :meth:`match` method can be used to compare two ``BIDSFile`` instances.
 
@@ -57,33 +60,45 @@ class BIDSFile(object):
         self.suffix      = suffix
 
 
-    def match(self, other):
+    def __str__(self):
+        """Return a strimg representation of this ``BIDSFile``. """
+        return 'BIDSFile({})'.format(self.filename)
+
+
+    def __repr__(self):
+        """Return a strimg representation of this ``BIDSFile``. """
+        return str(self)
+
+
+    def match(self, other, suffix=True):
         """Compare this ``BIDSFile`` to ``other``.
 
-        :arg other: ``BIDSFile`` to compare
-        :returns:   ``True`` if ``self.suffix == other.suffix`` and if
-                    all of the entities in ``other`` are present in ``self``,
-                    ``False`` otherwise.
+        :arg other:  ``BIDSFile`` to compare
+
+        :arg suffix: Defaults to ``True``. If ``False``, the comparison
+                     is made solely on the entity values.
+
+        :returns:    ``True`` if ``self.suffix == other.suffix`` (unless
+                     ``suffix`` is ``False``) and if all of the entities in
+                     ``other`` are present in ``self``, ``False`` otherwise.
         """
 
-        suffix   = self.suffix == other.suffix
+        suffix   = (not suffix) or (self.suffix == other.suffix)
         entities = True
 
         for key, value in other.entities.items():
-            entities = entities and self.entities.get(key, None) == value
+            entities = entities and (self.entities.get(key, None) == value)
 
         return suffix and entities
 
 
 def parseFilename(filename):
-    """Parses a BIDS-like file name. The file name must consist of zero or more
-    "entities" (alpha-numeric ``name-value`` pairs), a "suffix", all separated
-    by underscores, and a regular file extension. For example, the following
-    file::
+    """Parses a BIDS-like file name, returning the entities and suffix encoded
+    in the name. See the :func:`isBIDSFile` function for an explanation of
+    what is considered to be a valid BIDS file name.
 
-        sub-01_ses-01_task-stim_bold.nii.gz
-
-    has suffix ``bold``, and entities ``sub=01``, ``ses=01`` and ``task=stim``.
+    .. note:: This function assumes that no period (``.``) characters occur in
+              the body of a BIDS filename.
 
     :returns: A tuple containing:
                - A dict containing the entities
@@ -97,7 +112,7 @@ def parseFilename(filename):
     suffix   = None
     entities = []
     filename = op.basename(filename)
-    filename = fslpath.removeExt(filename, ['.nii', '.nii.gz', '.json'])
+    filename = fslpath.removeExt(filename, firstDot=True)
     parts    = filename.split('_')
 
     for part in parts[:-1]:
@@ -142,13 +157,26 @@ def inBIDSDir(filename):
 def isBIDSFile(filename, strict=True):
     """Returns ``True`` if ``filename`` looks like a BIDS image or JSON file.
 
+    A BIDS file name must consist of zero or more "entities" (alpha-numeric
+    ``name-value`` pairs), a "suffix", all separated by underscores, and a
+    regular file extension. For example, the following file::
+
+        sub-01_ses-01_task-stim_bold.nii.gz
+
+    has suffix ``bold``, entities ``sub=01``, ``ses=01`` and ``task=stim``, and
+    extension ``.nii.gz``.
+
     :arg filename: Name of file to check
     :arg strict:   If ``True`` (the default), the file must be within a BIDS
                    dataset directory, as defined by :func:`inBIDSDir`.
     """
 
+    # Zero or more entities because sidecar files
+    # do not necessarily need to contain any
+    # entities (e.g. ``T1w.json`` is valid)
+
     name    = op.basename(filename)
-    pattern = r'([a-z0-9]+-[a-z0-9]+_)*([a-z0-9])+\.(nii|nii\.gz|json)'
+    pattern = r'([a-z0-9]+-[a-z0-9]+_)*([a-z0-9])+\.(.+)'
     flags   = re.ASCII | re.IGNORECASE
     match   = re.fullmatch(pattern, name, flags) is not None
 
@@ -159,7 +187,7 @@ def isBIDSFile(filename, strict=True):
 def loadMetadataFile(filename):
     """Load ``filename`` (assumed to be JSON), returning its contents. """
     with open(filename, 'rt') as f:
-        return json.load(f)
+        return json.load(f, strict=False)
 
 
 def loadMetadata(filename):
@@ -170,7 +198,7 @@ def loadMetadata(filename):
                    ``filename``
     """
 
-    filename  = op.realpath(op.abspath(filename))
+    filename  = op.abspath(filename)
     bfile     = BIDSFile(filename)
     dirname   = op.dirname(filename)
     prevdir   = filename

fsl/utils/cache.py

@@ -16,10 +16,9 @@ class Expired(Exception):
     """``Exception`` raised by the :meth:`Cache.get` metho when an attempt is
     made to access a cache item that has expired.
     """
-    pass
 
 
-class CacheItem(object):
+class CacheItem:
     """Internal container class used to store :class:`Cache` items. """
 
     def __init__(self, key, value, expiry=0):
@@ -29,7 +28,7 @@ class CacheItem(object):
         self.storetime = time.time()
 
 
-class Cache(object):
+class Cache:
     """The ``Cache`` is a simple in-memory cache built on a
     ``collections.OrderedDict``. The ``Cache`` class has the following
     features:
@@ -150,6 +149,23 @@ class Cache(object):
         return key in self.__cache
 
 
+    def keys(self):
+        """Return all keys in the cache. """
+        return self.__cache.keys()
+
+
+    def values(self):
+        """Return all values in the cache. """
+        for item in self.__cache.values():
+            yield item.value
+
+
+    def items(self):
+        """Return all (key, value) pairs in the cache. """
+        for key, item in self.__cache.items():
+            yield key, item.value
+
+
     def __parseDefault(self, *args, **kwargs):
         """Used by the :meth:`get` method. Parses the ``default`` argument,
         which may be specified as either a positional or keyword argumnet.

fsl/utils/ensure.py

@@ -13,8 +13,7 @@ that some condition is met.
    ensureIsImage
 """
 
-
-import six
+import pathlib
 
 import nibabel as nib
 
@@ -24,7 +23,7 @@ import fsl.data.image as fslimage
 def ensureIsImage(img):
     """Ensures that the given ``img`` is an in-memory ``nibabel`` object.
     """
-    if isinstance(img, six.string_types):
+    if isinstance(img, (str, pathlib.Path)):
         img = fslimage.addExt(img)
         img = nib.load(img)
     return img