From f162bb5b987537acdf75997c65fcf2129052c4ae Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 5 Oct 2017 18:18:30 +0100
Subject: [PATCH] Image.resample uses scipy.ndimage.affine_transform instead of
 zoom, because the former aligns on the centre of voxel (0, 0, 0), whereas the
 latter aligns on the corner.

---
 fsl/data/image.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/fsl/data/image.py b/fsl/data/image.py
index 0f8a9995d..c85d3ab61 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -1137,7 +1137,7 @@ class Image(Nifti):
                    - A ``numpy`` array of shape ``shape``, containing an
                      interpolated copy of the data in this ``Image``.
 
-                   - A ``numpy`` array of shape ``(4, 4``), containing the
+                   - A ``numpy`` array of shape ``(4, 4)``, containing the
                      adjusted voxel-to-world transformation for the resampled
                      data.
         """
@@ -1152,13 +1152,21 @@ class Image(Nifti):
         if tuple(data.shape) != tuple(shape):
 
             oldShape = data.shape
-            zooms    = [float(shape[i]) / data.shape[i] for i in range(ndims)]
-            data     = ndimage.zoom(data, zooms, **kwargs)
-
+            zooms    = [float(data.shape[i]) / shape[i] for i in range(ndims)]
+            zooms    = np.diag(zooms)
+
+            data = ndimage.affine_transform(data,
+                                            zooms,
+                                            output_shape=shape,
+                                            **kwargs)
+
+            # 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)]
-            offset   = [(s - 1) / 2.0  for s in scale]
-            scale    = transform.scaleOffsetXform(scale, offset)
+            scale    = transform.scaleOffsetXform(scale, 0)
             xform    = transform.concat(self.voxToWorldMat, scale)
         else:
             xform = self.voxToWorldMat
-- 
GitLab