diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 63c63cafd6b800500622fe1c18b293e8b5d6d9a1..019e6e6eff6f74ab7c951eeb9f32b435507fb3fb 100644
--- a/CHANGELOG.rst
+++ b/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)
-----------------------------------
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 4cde726519f32c255c17a37dfc27f9452f4f67da..39f1f7facdd789e197340177f8e73dea77507c31 100644
--- a/fsl/data/image.py
+++ b/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
diff --git a/fsl/transform/affine.py b/fsl/transform/affine.py
index e0364bf94030c76d6df99a5222ce6439bed8df38..c7fb71f2a0571f8a1318a7e7b84749cfa498a9f2 100644
--- a/fsl/transform/affine.py
+++ b/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
diff --git a/fsl/utils/image/resample.py b/fsl/utils/image/resample.py
index 50dc0c4022cc5b32f763d8afd43f3d7fb31fba47..7679c56b02a4b31ff91e4fd3b499f63015d38792 100644
--- a/fsl/utils/image/resample.py
+++ b/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, :]
diff --git a/tests/test_image.py b/tests/test_image.py
index ace4dcb510006bc71acf2c9aceba4493d2944c9c..3bb77d01a3d8f6e7d559fc79796db7097c04615e 100644
--- a/tests/test_image.py
+++ b/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))
diff --git a/tests/test_transform/test_affine.py b/tests/test_transform/test_affine.py
index 7a69b86e4d79531f53d52be991ba5afe994b248a..79a17fba6f9a3895b124b06f279be1c9b9f1b0f9 100644
--- a/tests/test_transform/test_affine.py
+++ b/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)))