From 56fa033fd378f71687d2f53fbef90c4082188f5a Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Fri, 6 Oct 2017 12:33:45 +0100
Subject: [PATCH] Resample method now has a smoothing option, because otherwise
interpolation has no effect when downsampling to a new shape where voxel
centres are aligned.
---
fsl/data/image.py | 71 +++++++++++++++++++++++++++++++++--------------
1 file changed, 50 insertions(+), 21 deletions(-)
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 7a1000540..f914e1094 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -1131,21 +1131,37 @@ class Image(Nifti):
self.notify(topic='saveState')
- def resample(self, shape, sliceobj=None, dtype=None, **kwargs):
+ def resample(self,
+ newShape,
+ sliceobj=None,
+ dtype=None,
+ order=1,
+ smooth=True):
"""Returns a copy of the data in this ``Image``, resampled to the
specified ``shape``.
- :arg shape: Desired shape
-
- :arg dtype: ``numpy`` data type of the resampled data. If ``None``,
- the :meth:`dtype` of this ``Image`` is used.
+ :arg newShape: Desired shape. May containg floating point values,
+ in which case the resampled image will have shape
+ ``round(newShape)``, but the voxel sizes will
+ have scales ``self.shape / newShape``.
:arg sliceobj: Slice into this ``Image``. If ``None``, the whole
image is resampled, and it is assumed that it has the
same number of dimensions as ``shape``.
- All other arguments are passed through to the ``scipy.ndimage.zoom``
- function.
+ :arg dtype: ``numpy`` data type of the resampled data. If ``None``,
+ the :meth:`dtype` of this ``Image`` 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 smooth: If ``True`` (the default), the data is smoothed before
+ being resampled, but only along axes which are being
+ down-sampled (i.e. where
+ ``newShape[i] < self.shape[i]``).
:returns: A tuple containing:
@@ -1160,29 +1176,42 @@ class Image(Nifti):
if sliceobj is None: sliceobj = slice(None)
if dtype is None: dtype = self.dtype
- ndims = len(shape)
- data = self[sliceobj]
- data = np.array(data, dtype=dtype, copy=False)
+ data = self[sliceobj]
+ data = np.array(data, dtype=dtype, copy=False)
+
+ oldShape = np.array(data.shape)
+ newShape = np.array(newShape)
+
+ if not np.all(np.isclose(oldShape, newShape)):
- if tuple(data.shape) != tuple(shape):
+ ratio = oldShape / newShape
+ newShape = np.array(np.round(newShape), dtype=np.int)
+ scale = transform.scaleOffsetXform(ratio, 0)
- oldShape = data.shape
- zooms = [float(data.shape[i]) / shape[i] for i in range(ndims)]
- zooms = np.diag(zooms)
+ # If interpolating and smoothing, we apply a
+ # gaussian filter along axes with a resampling
+ # ratio greater than 1.1. We do this so that
+ # interpolation has an effect when down-sampling
+ # to a resolution where the voxel centres are
+ # aligned (as otherwise any interpolation regime
+ # will be equivalent to nearest neighbour). This
+ # more-or-less mimics the behaviour of FLIRT.
+ if order > 0 and smooth:
+ sigma = np.array(ratio)
+ sigma[ratio < 1.1] = 0
+ sigma[ratio >= 1.1] *= 0.425
+ data = ndimage.gaussian_filter(data, sigma)
data = ndimage.affine_transform(data,
- zooms,
- output_shape=shape,
- **kwargs)
+ scale[:3, :3],
+ output_shape=newShape,
+ order=order)
# Construct an affine transform which
# puts the resampled image into the
# same world coordinate system as this
# image.
- newShape = data.shape
- scale = [os / float(ns) for os, ns in zip(oldShape, newShape)]
- scale = transform.scaleOffsetXform(scale, 0)
- xform = transform.concat(self.voxToWorldMat, scale)
+ xform = transform.concat(self.voxToWorldMat, scale)
else:
xform = self.voxToWorldMat
--
GitLab