Merge branch 'enh/nifti-adjust' into 'master'

Enh/nifti adjust

Closes #368

See merge request fsl/fslpy!182

CHANGELOG.rst

@@ -7,6 +7,16 @@ order.
 -------------------------
 
 
+Added
+^^^^^
+
+
+* New :meth:`.Nifti.adjust` method, for creating a copy of a :class:`.Nifti`
+  header with adjusted shape, pixdims, and affine. This can be useful for
+  creating a resampling reference.
+* New :func:`.affine.rescale` function, for adjusting a scaling matrix.
+
+
 Fixed
 ^^^^^
 
@@ -14,6 +24,14 @@ Fixed
 * Improved the algorithm used by the :func:`.mesh.needsFixing` function.
 
 
+Deprecated
+^^^^^^^^^^
+
+
+* :func:`.calculateMatrix` - its functionality has been moved to the
+  :func:`.affine.rescale` function.
+
+
 2.7.0 (Wednesday 6th November 2019)
 -----------------------------------
 

fsl/data/image.py

@@ -196,7 +196,7 @@ class Nifti(notifier.Notifier, meta.Meta):
 
     A ``Nifti`` instance expects to be passed either a
     ``nibabel.nifti1.Nifti1Header`` or a ``nibabel.nifti2.Nifti2Header``, but
-    can als encapsulate a ``nibabel.analyze.AnalyzeHeader``. In this case:
+    can also encapsulate a ``nibabel.analyze.AnalyzeHeader``. In this case:
 
       - The image voxel orientation is assumed to be R->L, P->A, I->S.
 
@@ -839,6 +839,73 @@ class Nifti(notifier.Notifier, meta.Meta):
         return code
 
 
+    def adjust(self, pixdim=None, shape=None, origin=None):
+        """Return a new ``Nifti`` object with the specified ``pixdim`` or
+        ``shape``. The affine of the new ``Nifti`` is adjusted accordingly.
+
+        Only one of ``pixdim`` or ``shape`` can be specified.
+
+        See :func:`.affine.rescale` for the meaning of the ``origin`` argument.
+
+        Only the spatial dimensions may be adjusted - use the functions in
+        the :mod:`.image.resample` module if you need to adjust non-spatial
+        dimensions.
+
+        :arg pixdim: New voxel dimensions
+        :arg shape:  New image shape
+        :arg origin: Voxel grid alignment - either ``'centre'`` (the default)
+                     or ``'corner'``
+        :returns:    A new ``Nifti`` object based on this one, with adjusted
+                     pixdims, shape and affine.
+        """
+
+        if ((pixdim is not None) and (shape is not None)) or \
+           ((pixdim is     None) and (shape is     None)):
+            raise ValueError('Exactly one of pixdim or '
+                             'shape must be specified')
+
+        if shape is not None: ndim = len(shape)
+        else:                 ndim = len(pixdim)
+
+        # We only allow adjustment of
+        # the spatial dimensions
+        if ndim != 3:
+            raise ValueError('Three dimensions must be specified')
+
+        oldShape  = np.array(self.shape[ :ndim])
+        oldPixdim = np.array(self.pixdim[:ndim])
+        newShape  = shape
+        newPixdim = pixdim
+
+        # if pixdims were specified,
+        # convert them into a shape,
+        # and vice versa
+        if newPixdim is not None:
+            newShape = oldShape * (oldPixdim / newPixdim)
+        else:
+            newPixdim = oldPixdim * (oldShape / newShape)
+
+        # Rescale the voxel-to-world affine
+        xform = affine.rescale(oldShape, newShape, origin)
+        xform = affine.concat(self.getAffine('voxel', 'world'), xform)
+
+        # Now that we've got our spatial
+        # scaling/offset matrix, pad the
+        # new shape/pixdims with those
+        # from any non-spatial dimensions
+        newShape  = list(newShape)  + list(self.shape[ 3:])
+        newPixdim = list(newPixdim) + list(self.pixdim[3:])
+
+        # And create the new header
+        # and we're away
+        header = self.header.copy()
+        header.set_data_shape(newShape)
+        header.set_zooms(     newPixdim)
+        header.set_sform(     xform)
+        header.set_qform(     xform)
+        return Nifti(header)
+
+
 class Image(Nifti):
     """Class which represents a NIFTI image. Internally, the image is
     loaded/stored using a :mod:`nibabel.nifti1.Nifti1Image` or

fsl/transform/affine.py

@@ -21,6 +21,7 @@ transformations. The following functions are available:
    axisAnglesToRotMat
    axisBounds
    rmsdev
+   rescale
 
 And a few more functions are provided for working with vectors:
 
@@ -582,3 +583,59 @@ def rmsdev(T1, T2, R=None, xc=None):
     erms = np.sqrt(erms)
 
     return erms
+
+
+def rescale(oldShape, newShape, origin=None):
+    """Calculates an affine matrix to use for resampling.
+
+    This function generates an affine transformation matrix that can be used
+    to resample an N-D array from ``oldShape`` to ``newShape`` using, for
+    example, ``scipy.ndimage.affine_transform``.
+
+    The matrix will contain scaling factors derived from the ``oldShape /
+    newShape`` ratio, and an offset determined by the ``origin``.
+
+    The default value for ``origin`` (``'centre'``) causes the corner voxel of
+    the output to have the same centre as the corner voxel of the input. If
+    the origin is ``'corner'``, we apply an offset which effectively causes
+    the voxel grid corners of the input and output to be aligned.
+
+    :arg oldShape: Shape of input data
+    :arg newShape: Shape to resample data to
+    :arg origin:   Voxel grid alignment - either ``'centre'`` (the default) or
+                   ``'corner'``
+    :returns:      An affine resampling matrix
+    """
+
+    if origin is None:
+        origin = 'centre'
+
+    oldShape = np.array(oldShape, dtype=np.float)
+    newShape = np.array(newShape, dtype=np.float)
+    ndim     = len(oldShape)
+
+    if len(oldShape) != len(newShape):
+        raise ValueError('Shape mismatch')
+
+    # shapes are the same - no rescaling needed
+    if np.all(np.isclose(oldShape, newShape)):
+        return np.eye(ndim + 1)
+
+    # Otherwise we calculate a scaling
+    # matrix from the old/new shape
+    # ratio, and specify an offset
+    # according to the origin
+    ratio = oldShape / newShape
+    scale = np.diag(ratio)
+
+    # Calculate an offset from the origin
+    if   origin == 'centre': offset = [0] * ndim
+    elif origin == 'corner': offset = (ratio - 1) / 2
+
+    # combine the scales and translations
+    # to form thte final affine
+    xform               = np.eye(ndim + 1)
+    xform[:ndim, :ndim] = scale
+    xform[:ndim, -1]    = offset
+
+    return xform

fsl/utils/image/resample.py

@@ -10,17 +10,15 @@ to resample an :class:`.Image` object to a different resolution.
 The :func:`resampleToPixdims` and :func:`resampleToReference` functions
 are convenience wrappers around :func:`resample`.
 
-The :func:`applySmoothing` and :func:`calculateMatrix` functions are
-sub-functions of :func:`resample`.
+The :func:`applySmoothing` function is a sub-function of :func:`resample`.
 """
 
 
-import collections.abc      as abc
-
 import numpy                as np
 import scipy.ndimage        as ndimage
 
 import fsl.transform.affine as affine
+import fsl.utils.deprecated as deprecated
 
 
 def resampleToPixdims(image, newPixdims, **kwargs):
@@ -204,12 +202,11 @@ def resample(image,
     # old/new shape ratio and the origin
     # setting.
     if matrix is None:
-        matrix = calculateMatrix(data.shape, newShape, origin)
+        matrix = affine.rescale(data.shape, newShape, origin)
 
-    # calculateMatrix will return None
-    # if it decides that the image
+    # identity matrix? the image
     # doesn't need to be resampled
-    if matrix is None:
+    if np.all(np.isclose(matrix, np.eye(len(newShape) + 1))):
         return data, image.voxToWorldMat
 
     newShape = np.array(np.round(newShape), dtype=np.int)
@@ -230,9 +227,9 @@ def resample(image,
     # Construct an affine transform which
     # puts the resampled image into the
     # same world coordinate system as this
-    # image. The calculateMatrix function
-    # might not return a 4x4 matrix, so we
-    # make sure it is valid.
+    # image. We may be working with >3D data,
+    # so here we discard the non-spatial
+    # parts of the resampling matrix
     if matrix.shape != (4, 4):
         rotmat         = matrix[:3, :3]
         offsets        = matrix[:3, -1]
@@ -286,53 +283,11 @@ def applySmoothing(data, matrix, newShape):
     return ndimage.gaussian_filter(data, sigma)
 
 
+@deprecated.deprecated('2.8.0', '3.0.0',
+                       'Use fsl.transform.affine.rescale instead')
 def calculateMatrix(oldShape, newShape, origin):
-    """Calculates an affine matrix to use for resampling.
-
-    Called by :func:`resample`.  The matrix will contain scaling factors
-    determined from the ``oldShape / newShape`` ratio, and an offset
-    determined from the ``origin``.
-
-    :arg oldShape: Shape of input data
-    :arg newShape: Shape to resample data to
-    :arg origin:   Voxel grid alignment - either ``'centre'`` or ``'corner'``
-    :returns:      An affine matrix that can be passed to
-                   ``scipy.ndimage.affine_transform``.
-    """
-
-    oldShape = np.array(oldShape, dtype=np.float)
-    newShape = np.array(newShape, dtype=np.float)
-
-    if np.all(np.isclose(oldShape, newShape)):
+    """Deprecated - use :func:`.affine.rescale` instead. """
+    xform = affine.rescale(oldShape, newShape, origin)
+    if np.all(np.isclose(xform, np.eye(len(oldShape) + 1))):
         return None
-
-    # Otherwise we calculate a
-    # scaling matrix from the
-    # old/new shape ratio, and
-    # specify an offset
-    # according to the origin
-    else:
-        ratio = oldShape / newShape
-        scale = np.diag(ratio)
-
-        # Calculate an offset from the
-        # origin - the default behaviour
-        # (centre) causes the corner voxel
-        # of the output to have the same
-        # centre as the corner voxel of
-        # the input. If the origin is
-        # 'corner', we apply an offset
-        # which effectively causes the
-        # voxel grids of the input and
-        # output to be aligned.
-        if   origin == 'centre': offset = 0
-        elif origin == 'corner': offset = list((ratio - 1) / 2)
-
-        if not isinstance(offset, abc.Sequence):
-            offset = [offset] * len(newShape)
-
-        # ndimage.affine_transform will accept
-        # a matrix of shape (ndim, ndim + 1)
-        matrix = np.hstack((scale, np.atleast_2d(offset).T))
-
-    return matrix
+    return xform[:-1, :]

tests/test_image.py

@@ -1430,3 +1430,28 @@ def test_loadMetadata():
         gotmeta = fslimage.loadMetadata(img)
 
         assert gotmeta == meta
+
+
+def test_adjust():
+
+    with tempdir():
+        make_image('image.nii', dims=(10, 10, 10), pixdims=(1, 1, 1))
+
+        img = fslimage.Image('image.nii')
+
+        assert img.sameSpace(img.adjust(pixdim=(1, 1, 1)))
+        assert img.sameSpace(img.adjust(shape=(10, 10, 10)))
+
+        adj = img.adjust(shape=(20, 20, 20))
+        assert adj.shape  == (20,  20,  20)
+        assert adj.pixdim == (0.5, 0.5, 0.5)
+
+        adj = img.adjust(pixdim=(0.5, 0.5, 0.5))
+        assert adj.shape  == (20,  20,  20)
+        assert adj.pixdim == (0.5, 0.5, 0.5)
+
+
+        adj  = img.adjust(shape=(8, 7, 11), origin='corner')
+        imgb = affine.axisBounds(img.shape, img.voxToWorldMat)
+        adjb = affine.axisBounds(adj.shape, adj.voxToWorldMat)
+        assert np.all(np.isclose(imgb, adjb, rtol=1e-5, atol=1e-5))

tests/test_transform/test_affine.py

@@ -528,3 +528,48 @@ def test_rmsdev():
         assert result < lastdist
 
         lastdist = result
+
+
+def test_rescale():
+
+    with pytest.raises(ValueError):
+        affine.rescale((5, 5), (10, 10, 10))
+
+    assert np.all(affine.rescale((5, 5),       (5, 5))       == np.eye(3))
+    assert np.all(affine.rescale((5, 5, 5),    (5, 5, 5))    == np.eye(4))
+    assert np.all(affine.rescale((5, 5, 5, 5), (5, 5, 5, 5)) == np.eye(5))
+
+    # (old shape, new shape, origin, expect)
+    tests = [
+        ((5, 5), (10, 10), 'centre', np.array([[0.5, 0,    0],
+                                               [0,   0.5,  0],
+                                               [0,   0,    1]])),
+        ((5, 5), (10, 10), 'corner', np.array([[0.5, 0,   -0.25],
+                                               [0,   0.5, -0.25],
+                                               [0,   0,    1]])),
+        ((5, 5, 5), (10, 10, 10), 'centre', np.array([[0.5, 0,    0,   0],
+                                                      [0,   0.5,  0,   0],
+                                                      [0,   0,    0.5, 0],
+                                                      [0,   0,    0,   1]])),
+        ((5, 5, 5), (10, 10, 10), 'corner', np.array([[0.5, 0,    0,   -0.25],
+                                                      [0,   0.5,  0,   -0.25],
+                                                      [0,   0,    0.5, -0.25],
+                                                      [0,   0,    0,    1]])),
+        ((5, 5, 5, 5), (10, 10, 10, 10), 'centre', np.array([[0.5, 0,    0,   0,   0],
+                                                             [0,   0.5,  0,   0,   0],
+                                                             [0,   0,    0.5, 0,   0],
+                                                             [0,   0,    0,   0.5, 0],
+                                                             [0,   0,    0,   0,   1]])),
+        ((5, 5, 5, 5), (10, 10, 10, 10), 'corner', np.array([[0.5, 0,    0,   0,   -0.25],
+                                                             [0,   0.5,  0,   0,   -0.25],
+                                                             [0,   0,    0.5, 0,   -0.25],
+                                                             [0,   0,    0,   0.5, -0.25],
+                                                             [0,   0,    0,   0,   1]])),
+    ]
+
+    for oldshape, newshape, origin, expect in tests:
+
+        got = affine.rescale(oldshape, newshape, origin)
+        assert np.all(np.isclose(got, expect))
+        got = affine.rescale(newshape, oldshape, origin)
+        assert np.all(np.isclose(got, affine.invert(expect)))