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

@@ -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