diff --git a/.ci/test_template.sh b/.ci/test_template.sh
index 00bc034b1b76c8db636817e86f3f7820150a250c..c3ca00d519a70e21135f6a2720ea94b328f563ad 100644
--- a/.ci/test_template.sh
+++ b/.ci/test_template.sh
@@ -48,11 +48,11 @@ xvfb-run python setup.py test --addopts="$TEST_OPTS tests/test_platform.py"
 # unintuitively, includes nobody)
 chmod -R a+w `pwd`
 cmd="source /test.venv/bin/activate && python setup.py test"
-cmd="$cmd --addopts='$TEST_OPTS tests/test_immv_imcp.py'"
+cmd="$cmd --addopts='$TEST_OPTS tests/test_scripts/test_immv_imcp.py tests/test_immv_imcp.py'"
 su -s /bin/bash -c "$cmd" nobody
 
 # All other tests can be run as normal.
-python setup.py test --addopts="$TEST_OPTS -m 'not longtest' --ignore=tests/test_idle.py --ignore=tests/test_platform.py --ignore=tests/test_immv_imcp.py"
+python setup.py test --addopts="$TEST_OPTS -m 'not longtest' --ignore=tests/test_idle.py --ignore=tests/test_platform.py --ignore=tests/test_immv_imcp.py --ignore=tests/test_scripts/test_immv_imcp.py"
 
 # Long tests are only run on release branches
 if [[ $CI_COMMIT_REF_NAME == v* ]]; then
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3c259ef294d3e2c5e9891356cf22d8ac34fc38f6..1088480c33d93f9196a17c4fadfa2c7735108f11 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,6 +2,49 @@ This document contains the ``fslpy`` release history in reverse chronological
 order.
 
 
+2.2.0 (Wednesday May 8th 2019)
+------------------------------
+
+
+Added
+^^^^^
+
+
+* New :mod:`.resample_image` script.
+* New :mod:`.resample` module (replacing the :func:`.Image.resample` method),
+  containing functions to resample an :class:`.Image`.
+* New :func:`.resample.resampleToPixdim` and
+  :func:`.resample.resampleToReference` functions, convenience wrappers around
+  :func:`.resample.resample`.
+* New :func:`.idle.block` function.
+
+
+Changed
+^^^^^^^
+
+
+* The :func:`.resample` function (formerly :meth:`.Image.resample`) now
+  accepts ``origin`` and ``matrix`` parameters, which can be used to adjust
+  the alignment of the voxel grids of the input and output images.
+* The :func:`.transform.decompose` function now accepts both ``(3, 3)``
+  and ``(4, 4)`` matrices.
+
+
+Fixed
+^^^^^
+
+
+* Minor fixes to some :mod:`.filetree` tree definitions.
+
+
+Deprecated
+^^^^^^^^^^
+
+
+* The :meth:`.Image.resample` method has been deprecated in favour of the
+  :func:`.resample.resample` function.
+
+
 2.1.0 (Saturday April 13th 2019)
 --------------------------------
 
diff --git a/fsl/data/atlases.py b/fsl/data/atlases.py
index 259b91fdf41bc84490a6ecf28344d09ce3ce8c19..6ec39a89701554f55c4a943cdd3315e02cf2b004 100644
--- a/fsl/data/atlases.py
+++ b/fsl/data/atlases.py
@@ -53,6 +53,7 @@ import numpy                              as np
 import fsl.data.image                     as fslimage
 import fsl.data.constants                 as constants
 from   fsl.utils.platform import platform as platform
+import fsl.utils.image.resample           as resample
 import fsl.utils.transform                as transform
 import fsl.utils.notifier                 as notifier
 import fsl.utils.settings                 as fslsettings
@@ -695,9 +696,8 @@ class Atlas(fslimage.Image):
         # for resampling, as it is most likely
         # that the mask is binary.
         try:
-            mask, xform = mask.resample(self.shape[:3],
-                                        dtype=np.float32,
-                                        order=0)
+            mask, xform = resample.resample(
+                mask, self.shape[:3], dtype=np.float32, order=0)
 
         except ValueError:
             raise MaskError('Mask has wrong number of dimensions')
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 24e129c8bfc1ecc42d56d5c77d241e3912c5bd67..651992c4dbc3401a5af75ee3e4b59ea2f153fc8e 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -41,7 +41,6 @@ import                      warnings
 
 import                      six
 import numpy             as np
-import scipy.ndimage     as ndimage
 
 import nibabel           as nib
 import nibabel.fileslice as fileslice
@@ -1166,96 +1165,12 @@ class Image(Nifti):
         self.notify(topic='saveState')
 
 
-    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 ``newShape``.
-
-        :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  ``newShape``. A
-                       :exc:`ValueError` is raised if this is not the case.
-
-        :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:
-
-                   - A ``numpy`` array of shape ``newShape``, containing
-                     an interpolated copy of the data in this ``Image``.
-
-                   - A ``numpy`` array of shape ``(4, 4)``, containing the
-                     adjusted voxel-to-world transformation for the spatial
-                     dimensions of the resampled data.
-        """
-
-        if sliceobj is None: sliceobj = slice(None)
-        if dtype    is None: dtype    = self.dtype
-
-        data = self[sliceobj]
-        data = np.array(data, dtype=dtype, copy=False)
-
-        oldShape = np.array(data.shape, dtype=np.float)
-        newShape = np.array(newShape,   dtype=np.float)
-
-        if len(oldShape) != len(newShape):
-            raise ValueError('Shapes don\'t match')
-
-        if not np.all(np.isclose(oldShape, newShape)):
-
-            ratio    = oldShape / newShape
-            newShape = np.array(np.round(newShape), dtype=np.int)
-            scale    = np.diag(ratio)
-
-            # 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,
-                                            scale,
-                                            output_shape=newShape,
-                                            order=order)
-
-            # Construct an affine transform which
-            # puts the resampled image into the
-            # same world coordinate system as this
-            # image.
-            scale = transform.scaleOffsetXform(ratio[:3], 0)
-            xform = transform.concat(self.voxToWorldMat, scale)
-        else:
-            xform = self.voxToWorldMat
-
-        return data, xform
+    @deprecated.deprecated('2.2.0', '3.0.0',
+                           'Use fsl.utils.image.resample instead.')
+    def resample(self, *args, **kwargs):
+        """Deprecated - use :func:`.image.resample` instead. """
+        from fsl.utils.image.resample import resample
+        return resample(self, *args, **kwargs)
 
 
     def __getitem__(self, sliceobj):
diff --git a/fsl/scripts/resample_image.py b/fsl/scripts/resample_image.py
new file mode 100644
index 0000000000000000000000000000000000000000..151c301bf3af108485a67f3d21fd65df9efedca2
--- /dev/null
+++ b/fsl/scripts/resample_image.py
@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+#
+# resample_image.py - Script to resample an image
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module defines the ``resample_image`` script, for resampling
+a NIfTI image.
+"""
+
+
+import textwrap as tw
+import             sys
+import             argparse
+
+import numpy    as np
+
+import fsl.utils.parse_data     as parse_data
+import fsl.utils.image.resample as resample
+import fsl.data.image           as fslimage
+
+
+ARGS = {
+    'input'     : ('input',),
+    'output'    : ('output',),
+    'shape'     : ('-s',  '--shape'),
+    'dim'       : ('-d',  '--dim'),
+    'reference' : ('-r',  '--reference'),
+    'interp'    : ('-i',  '--interp'),
+    'origin'    : ('-o',  '--origin'),
+    'dtype'     : ('-dt', '--dtype'),
+    'smooth'    : ('-n',  '--nosmooth')}
+
+
+OPTS = {
+    'input'     : dict(type=parse_data.Image),
+    'output'    : dict(type=parse_data.ImageOut),
+    'reference' : dict(type=parse_data.Image, metavar='IMAGE'),
+    'shape'     : dict(type=int,   nargs=3, metavar=('X', 'Y', 'Z')),
+    'dim'       : dict(type=float, nargs=3, metavar=('X', 'Y', 'Z')),
+    'interp'    : dict(choices=('nearest', 'linear', 'cubic'),
+                       default='linear'),
+    'origin'    : dict(choices=('centre', 'corner'), default='centre'),
+    'dtype'     : dict(choices=('char', 'short', 'int', 'float', 'double')),
+    'smooth'    : dict(dest='smooth', action='store_false')}
+
+
+HELPS = {
+    'input'     : 'Input image',
+    'output'    : 'Output image',
+    'shape'     : 'Output shape',
+    'dim'       : 'Output voxel dimensions',
+    'reference' : 'Resample input to the space of this reference image'
+                  '(overrides --origin)',
+    'interp'    : 'Interpolation (default: linear)',
+    'origin'    : 'Resampling origin (default: centre)',
+    'dtype'     : 'Data type (default: data type of input image)',
+    'smooth'    : 'Do not smooth image when downsampling'}
+
+
+DESC = tw.dedent("""
+Resample an image to different dimensions.
+""").strip()
+
+
+DEST_DESC = tw.dedent("""
+Specify the resampling destination space using one of the following
+options. Note that the --reference option will cause the field-of-view
+of the input image to be changed to that of the reference image.
+""").strip()
+
+
+USAGE = 'resample_image (--shape|--dim|--reference) [options] input output'
+
+
+INTERPS = {'nearest' : 0,
+           'linear'  : 1,
+           'cubic'   : 3}
+DTYPES  = {'char'    : np.uint8,
+           'short'   : np.int16,
+           'int'     : np.int32,
+           'float'   : np.float32,
+           'double'  : np.float64}
+
+
+def parseArgs(argv):
+    """Parses command-line arguments.
+
+    :arg argv: Sequence of command-line arguments
+    :returns:  An ``argparse.Namespace`` object containing parsed arguments.
+    """
+
+    parser = argparse.ArgumentParser(prog='resample_image',
+                                     usage=USAGE,
+                                     description=DESC)
+    dest   = parser.add_argument_group('Resampling destination', DEST_DESC)
+    dest   = dest.add_mutually_exclusive_group(required=True)
+
+    for a in ('input', 'output', 'interp', 'origin',
+              'dtype', 'smooth'):
+        parser.add_argument(*ARGS[a], help=HELPS[a], **OPTS[a])
+
+    for a in ('shape', 'dim', 'reference'):
+        dest.add_argument(*ARGS[a], help=HELPS[a], **OPTS[a])
+
+    if len(argv) == 0:
+        parser.print_help()
+        sys.exit(0)
+
+    args        = parser.parse_args(argv)
+    args.interp = INTERPS[   args.interp]
+    args.dtype  = DTYPES.get(args.dtype, args.input.dtype)
+
+    return args
+
+
+def main(argv=None):
+    """Entry point for ``resample_image``. Parses arguments, resamples the
+    input image, and saves it to the specified output file.
+
+    :arg argv: Sequence of command-line arguments. If not provided, taken
+               from ``sys.argv``.
+    """
+
+    if argv is None:
+        argv = sys.argv[1:]
+
+    args      = parseArgs(argv)
+    reskwargs = {
+        'dtype'  : args.dtype,
+        'order'  : args.interp,
+        'smooth' : args.smooth,
+        'origin' : args.origin}
+
+    # One of these is guaranteed to be set
+    if args.shape is not None:
+        func    = resample.resample
+        resargs = (args.input, args.shape)
+
+    elif args.dim is not None:
+        func    = resample.resampleToPixdims
+        resargs = (args.input, args.dim)
+
+    elif args.reference is not None:
+        func    = resample.resampleToReference
+        resargs = (args.input, args.reference)
+
+    resampled, xform = func(*resargs, **reskwargs)
+
+    if args.reference is None:
+        hdr = args.input.header
+    else:
+        hdr   = args.reference.header
+        xform = None
+
+    resampled = fslimage.Image(resampled, xform=xform, header=hdr)
+    resampled.save(args.output)
+
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())
diff --git a/fsl/utils/idle.py b/fsl/utils/idle.py
index cdd86a01614ce9435dcc8b0ad5673503a874140d..1425ba191dd6d954c074369088c0a33d532dfeff 100644
--- a/fsl/utils/idle.py
+++ b/fsl/utils/idle.py
@@ -19,6 +19,7 @@ Idle tasks
 .. autosummary::
    :nosignatures:
 
+   block
    idle
    idleWhen
    inIdle
@@ -377,6 +378,29 @@ def cancelIdle(taskName):
     _idleQueueDict[taskName].timeout = -1
 
 
+def block(secs, delta=0.01):
+    """Blocks for the specified number of seconds, yielding to the main ``wx``
+    loop.
+
+    If ``wx`` is not available, or a ``wx`` application is not running, this
+    function is equivalent to ``time.sleep(secs)``.
+
+    :arg secs:  Time in seconds to block
+    :arg delta: Time in seconds to sleep between successive yields to ``wx``.
+    """
+
+    from fsl.utils.platform import platform as fslplatform
+
+    if not fslplatform.haveGui:
+        time.sleep(secs)
+    else:
+        import wx
+        start = time.time()
+        while (time.time() - start) < secs:
+            wx.YieldIfNeeded()
+            time.sleep(delta)
+
+
 def idle(task, *args, **kwargs):
     """Run the given task on a ``wx.EVT_IDLE`` event.
 
@@ -556,7 +580,7 @@ def idleWhen(func, condition, *args, **kwargs):
 
 def wait(threads, task, *args, **kwargs):
     """Creates and starts a new ``Thread`` which waits for all of the ``Thread``
-    instances to finsih (by ``join``ing them), and then runs the given
+    instances to finish (by ``join``ing them), and then runs the given
     ``task`` via :func:`idle`.
 
     If the ``direct`` parameter is ``True``, or a ``wx.App`` is not running,
diff --git a/fsl/utils/image/__init__.py b/fsl/utils/image/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..8e22533648b6ff1c4125b31e0375b594f937164b
--- /dev/null
+++ b/fsl/utils/image/__init__.py
@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+#
+# __init__.py - The fsl.utils.image package
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""The :mod:`fsl.utils.image` oackage contains algorithms and utilities for
+manipulating and working with :class:`.Image` objects.
+
+The following modules are available:
+
+.. autosumary::
+   :nosignature
+
+   .image.resample
+"""
diff --git a/fsl/utils/image/resample.py b/fsl/utils/image/resample.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc49a621dca6a5759421e8353578fbb7038eca8e
--- /dev/null
+++ b/fsl/utils/image/resample.py
@@ -0,0 +1,273 @@
+#!/usr/bin/env python
+#
+# resample.py - The resample functino
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module defines the :func:`resample` function, which can be used
+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`.
+"""
+
+
+import collections.abc     as abc
+
+import numpy               as np
+import scipy.ndimage       as ndimage
+
+import fsl.utils.transform as transform
+
+
+def resampleToPixdims(image, newPixdims, **kwargs):
+    """Resample ``image`` so that it has the specified voxel dimensions.
+
+    This is a wrapper around :func:`resample` - refer to its documenttion
+    for details on the other arguments and the return values.
+
+    :arg image:   :class:`.Image` to resample
+    :arg pixdims: New voxel dimensions to resample ``image`` to.
+    """
+    newPixdims = np.array(newPixdims)
+    oldShape   = np.array(image.shape)
+    oldPixdims = np.array(image.pixdim)
+    newShape   = oldShape * (oldPixdims / newPixdims)
+    return resample(image, newShape, **kwargs)
+
+
+def resampleToReference(image, reference, **kwargs):
+    """Resample ``image`` into the space of the ``reference``.
+
+    This is a wrapper around :func:`resample` - refer to its documenttion
+    for details on the other arguments and the return values.
+
+    :arg image:     :class:`.Image` to resample
+    :arg reference: :class:`.Nifti` defining the space to resample ``image``
+                    into
+    """
+
+    kwargs['mode']     = kwargs.get('mode', 'constant')
+    kwargs['newShape'] = reference.shape
+    kwargs['matrix']   = transform.concat(image.worldToVoxMat,
+                                          reference.voxToWorldMat)
+    return resample(image, **kwargs)
+
+
+def resample(image,
+             newShape,
+             sliceobj=None,
+             dtype=None,
+             order=1,
+             smooth=True,
+             origin='centre',
+             matrix=None,
+             mode='nearest',
+             cval=0):
+    """Returns a copy of the data in the ``image``, resampled to the specified
+    ``newShape``.
+
+    The space that the image is resampled into can be defined in one of the
+    following ways, in decreasing order of precedence:
+
+      1. If a ``matrix`` is provided, it is applied to the voxel coordinates
+         when retrieving values from the ``image``
+
+      2. Otherwise the image is simply scaled according to the ratio calculated
+         by ``image.shape / newShape``. In this case the ``origin`` argument
+         may be used to adjust the alignemnt of the original and resampled
+         voxel grids.
+
+    See the ``scipy.ndimage.affine_transform`` function for more details,
+    particularly on the ``order``, ``matrix``, ``mode`` and
+    ``cval`` arguments.
+
+    :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`` (unless ``matrix`` is specified).
+
+    :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  ``newShape``. A
+                   :exc:`ValueError` is raised if this is not the case.
+
+    :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]``).
+
+    :arg origin:   ``'centre'`` (the default) or ``'corner'``. ``'centre'``
+                   resamples the image such that the centre of the corner
+                   voxels of this image and the resampled data are
+                   aligned. ``'corner'`` resamples the image such that
+                   the corner of the corner voxels are aligned (and
+                   therefore the voxel grids are aligned).
+                   Ignored if ``offset`` or ``matrix`` is specified.
+
+    :arg matrix:   Arbitrary affine transformation matrix to apply to the
+                   voxel coordinates of ``image`` when resampling.
+
+    :arg mode:     How to handle regions which are outside of the image FOV.
+                   Defaults to `''nearest'``.
+
+    :arg cval:     Constant value to use when ``mode='constant'``.
+
+    :returns: A tuple containing:
+
+               - A ``numpy`` array of shape ``newShape``, containing
+                 an interpolated copy of the data in this ``Image``.
+
+               - A ``numpy`` array of shape ``(4, 4)``, containing the
+                 adjusted voxel-to-world transformation for the spatial
+                 dimensions of the resampled data.
+    """
+
+    if sliceobj is None:     sliceobj = slice(None)
+    if dtype    is None:     dtype    = image.dtype
+    if origin   == 'center': origin   = 'centre'
+
+    if origin not in ('centre', 'corner'):
+        raise ValueError('Invalid value for origin: {}'.format(origin))
+
+    data = np.array(image[sliceobj], dtype=dtype, copy=False)
+
+    if len(data.shape) != len(newShape):
+        raise ValueError('Data dimensions do not match new shape: '
+                         'len({}) != len({})'.format(data.shape, newShape))
+
+    # If matrix not provided, calculate
+    # a scaling/offset matrix from the
+    # old/new shape ratio and the origin
+    # setting.
+    if matrix is None:
+        matrix = calculateMatrix(data.shape, newShape, origin)
+
+    # calculateMatrix will return None
+    # if it decides that the image
+    # doesn't need to be resampled
+    if matrix is None:
+        return data, image.voxToWorldMat
+
+    newShape = np.array(np.round(newShape), dtype=np.int)
+
+    # Apply smoothing if requested,
+    # and if not using nn interp
+    if order > 0 and smooth:
+        data = applySmoothing(data, matrix, newShape)
+
+    # Do the resample thing
+    data = ndimage.affine_transform(data,
+                                    matrix,
+                                    output_shape=newShape,
+                                    order=order,
+                                    mode=mode,
+                                    cval=cval)
+
+    # 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.
+    if matrix.shape != (4, 4):
+        matrix = np.vstack((matrix[:3, :4], [0, 0, 0, 1]))
+    matrix = transform.concat(image.voxToWorldMat, matrix)
+
+    return data, matrix
+
+
+def applySmoothing(data, matrix, newShape):
+    """Called by the :func:`resample` function.
+
+    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.
+
+    See the ``scipy.ndimage.gaussian_filter`` function for more details.
+
+    :arg data:     Data to be smoothed.
+    :arg matrix:   Affine matrix to be used during resampling. The voxel
+                   scaling factors are extracted from this.
+    :arg newShape: Shape the data is to be resampled into.
+    :returns:      A smoothed copy of ``data``.
+    """
+
+    ratio = transform.decompose(matrix[:3, :3])[0]
+
+    if len(newShape) > 3:
+        ratio = np.concatenate((
+            ratio,
+            [float(o) / float(s)
+             for o, s in zip(data.shape[3:], newShape[3:])]))
+
+    sigma                = np.array(ratio)
+    sigma[ratio <  1.1]  = 0
+    sigma[ratio >= 1.1] *= 0.425
+
+    return ndimage.gaussian_filter(data, sigma)
+
+
+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)):
+        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
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index f574eaa764a8e337cb7616f7eaef254e0dfd3a9d..8c9be2c440a139851371d9f1e86454915b11fb7a 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -178,7 +178,7 @@ def decompose(xform, angles=True):
     It is assumed that the given transform has no perspective components. Any
     shears in the affine are discarded.
 
-    :arg xform:  A ``(4, 4)`` affine transformation matrix.
+    :arg xform:  A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
 
     :arg angles: If ``True`` (the default), the rotations are returned
                  as axis-angles, in radians. Otherwise, the rotation matrix
@@ -187,7 +187,8 @@ def decompose(xform, angles=True):
     :returns: The following:
 
                - A sequence of three scales
-               - A sequence of three translations
+               - A sequence of three translations (all ``0`` if ``xform``
+                 was a ``(3, 3)`` matrix)
                - A sequence of three rotations, in radians. Or, if
                  ``angles is False``, a rotation matrix.
     """
@@ -198,9 +199,13 @@ def decompose(xform, angles=True):
     # The next step is to extract the translations. This is trivial;
     # we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
     # point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
-    xform        = xform.T
-    translations = xform[ 3, :3]
-    xform        = xform[:3, :3]
+    xform = xform.T
+
+    if xform.shape == (4, 4):
+        translations = xform[ 3, :3]
+        xform        = xform[:3, :3]
+    else:
+        translations = np.array([0, 0, 0])
 
     M1 = xform[0]
     M2 = xform[1]
@@ -261,7 +266,7 @@ def decompose(xform, angles=True):
     if angles: rotations = rotMatToAxisAngles(R.T)
     else:      rotations = R.T
 
-    return [sx, sy, sz], translations, rotations
+    return np.array([sx, sy, sz]), translations, rotations
 
 
 def rotMatToAffine(rotmat, origin=None):
diff --git a/setup.py b/setup.py
index 1383c0e1228c3093d0de4433853e5f1adf2ec26d..0873bffe228abb22dc7a04647fe6a7b4ff93b5bc 100644
--- a/setup.py
+++ b/setup.py
@@ -119,12 +119,13 @@ setup(
 
     entry_points={
         'console_scripts' : [
-            'immv          = fsl.scripts.immv:main',
-            'imcp          = fsl.scripts.imcp:main',
-            'imglob        = fsl.scripts.imglob:main',
-            'atlasq        = fsl.scripts.atlasq:main',
-            'atlasquery    = fsl.scripts.atlasq:atlasquery_emulation',
-            'fsl_ents = fsl.scripts.fsl_ents:main',
+            'immv           = fsl.scripts.immv:main',
+            'imcp           = fsl.scripts.imcp:main',
+            'imglob         = fsl.scripts.imglob:main',
+            'atlasq         = fsl.scripts.atlasq:main',
+            'atlasquery     = fsl.scripts.atlasq:atlasquery_emulation',
+            'fsl_ents       = fsl.scripts.fsl_ents:main',
+            'resample_image = fsl.scripts.resample_image:main',
         ]
     }
 )
diff --git a/tests/test_atlases.py b/tests/test_atlases.py
index 80be87c74736f39cfa4134a3efa4e804d2ce0db3..c2d257f601da2eaf61b5de4eb6f4cd3e148d8d7d 100644
--- a/tests/test_atlases.py
+++ b/tests/test_atlases.py
@@ -17,9 +17,10 @@ import mock
 import pytest
 
 import tests
-import fsl.utils.transform as transform
-import fsl.data.atlases    as atlases
-import fsl.data.image      as fslimage
+import fsl.utils.transform      as transform
+import fsl.utils.image.resample as resample
+import fsl.data.atlases         as atlases
+import fsl.data.image           as fslimage
 
 
 datadir = op.join(op.dirname(__file__), 'testdata')
@@ -286,7 +287,7 @@ def test_prepareMask():
             np.array(np.random.random(ashape), dtype=np.float32),
             xform=atlas.voxToWorldMat)
 
-        goodmask2, xf = goodmask1.resample(m2shape)
+        goodmask2, xf = resample.resample(goodmask1, m2shape)
         goodmask2     = fslimage.Image(goodmask2, xform=xf)
 
         wrongdims     = fslimage.Image(
diff --git a/tests/test_atlases_query.py b/tests/test_atlases_query.py
index f90c8d40618b78ea7f33e0c8168623fcdac31d7c..693366a2e65000685196f4c7067304cd8a59fad9 100644
--- a/tests/test_atlases_query.py
+++ b/tests/test_atlases_query.py
@@ -14,6 +14,7 @@ import                    pytest
 import fsl.data.atlases    as fslatlases
 import fsl.data.image      as fslimage
 import fsl.utils.transform as transform
+import fsl.utils.image.resample as resample
 import fsl.utils.cache     as cache
 
 from . import (testdir, make_random_mask)
@@ -256,7 +257,8 @@ def _gen_mask_query(atlas, qtype, qin, maskres):
         # aggresively to make sure there
         # is no overlap between the different
         # resolutions
-        mask, xform = mask.resample(a.shape[:3], dtype=np.float32, order=1)
+        mask, xform = resample.resample(
+            mask, a.shape[:3], dtype=np.float32, order=1)
 
         mask[mask   < 1.0] = 0
         mask[a_zmask == 0] = 0
@@ -279,7 +281,8 @@ def _eval_mask_query(atlas, query, qtype, qin):
     if maskres == res:
         rmask = mask[:]
     else:
-        rmask = mask.resample(atlas.shape[:3], dtype=np.float32, order=0)[0]
+        rmask = resample.resample(
+            mask, atlas.shape[:3], dtype=np.float32, order=0)[0]
 
     rmask = np.array(rmask, dtype=np.bool)
 
diff --git a/tests/test_idle.py b/tests/test_idle.py
index 13a8e3c98d679fb1fe5dcb10044fe163fe6f852b..430e49fa06522681bae391320e898fb12180cbdb 100644
--- a/tests/test_idle.py
+++ b/tests/test_idle.py
@@ -157,6 +157,29 @@ def _test_idleTimeout():
     assert idle.getIdleTimeout() == default
 
 
+@pytest.mark.wxtest
+def test_block_with_gui():    _run_with_wx(   _test_block)
+def test_block_without_gui(): _run_without_wx(_test_block)
+def _test_block():
+
+    called = [False]
+
+    if fslplatform.haveGui:
+        import wx
+        def idlefunc():
+            called[0] = True
+        wx.CallLater(1000, idlefunc)
+
+    start = time.time()
+
+    idle.block(2)
+    end = time.time()
+    assert (end - start) >= 2
+
+    if fslplatform.haveGui:
+        assert called[0]
+
+
 @pytest.mark.wxtest
 def test_idle():
 
diff --git a/tests/test_image.py b/tests/test_image.py
index 3290ed3619e6d13c713977cf851dd80bd37843c8..85dde7c66d44dd76bc137856f78de72163df2d39 100644
--- a/tests/test_image.py
+++ b/tests/test_image.py
@@ -1065,127 +1065,6 @@ def _test_Image_save(imgtype):
             img2 = None
 
 
-def test_image_resample(seed):
-
-    with tempdir() as td:
-
-        fname = op.join(td, 'test.nii')
-
-        # Random base image shapes
-        for i in range(25):
-
-            shape = np.random.randint(5, 50, 3)
-            make_random_image(fname, shape)
-            img = fslimage.Image(fname, mmap=False)
-
-            # bad shape
-            with pytest.raises(ValueError):
-                img.resample((10, 10))
-            with pytest.raises(ValueError):
-                img.resample((10, 10, 10, 10))
-
-            # resampling to the same shape should be a no-op
-            samei, samex = img.resample(shape)
-            assert np.all(samei == img[:])
-            assert np.all(samex == img.voxToWorldMat)
-
-            # Random resampled image shapes
-            for j in range(10):
-
-                rshape        = np.random.randint(5, 50, 3)
-                resampled, xf = img.resample(rshape, order=0)
-
-                img.save('base.nii.gz')
-                fslimage.Image(resampled, xform=xf,
-                               mmap=False).save('res.nii.gz')
-
-                assert tuple(resampled.shape) == tuple(rshape)
-
-                # We used nearest neighbour interp, so the
-                # values in the resampled image should match
-                # corresponding values in the original. Let's
-                # check some whynot.
-                restestcoords = np.array([
-                    np.random.randint(0, rshape[0], 100),
-                    np.random.randint(0, rshape[1], 100),
-                    np.random.randint(0, rshape[2], 100)]).T
-
-                resx,  resy,  resz  = restestcoords.T
-                resvals  = resampled[resx, resy, resz]
-
-                res2orig = transform.concat(img.worldToVoxMat, xf)
-
-                origtestcoords = transform.transform(restestcoords, res2orig)
-
-                # remove any coordinates which are out of
-                # bounds in the original image space, or
-                # are right on a voxel boundary (where the
-                # nn interp could have gone either way), or
-                # have value == 0 in the resampled space.
-                out = ((origtestcoords < 0)            |
-                       (origtestcoords >= shape - 0.5) |
-                       (np.isclose(np.modf(origtestcoords)[0], 0.5)))
-                out = np.any(out, axis=1) | (resvals == 0)
-
-                origtestcoords = np.array(origtestcoords.round(), dtype=np.int)
-
-                origtestcoords = origtestcoords[~out, :]
-                restestcoords  = restestcoords[ ~out, :]
-
-                resx,  resy,  resz  = restestcoords.T
-                origx, origy, origz = origtestcoords.T
-
-                origvals = img[:][origx, origy, origz]
-                resvals  = resampled[resx, resy, resz]
-
-                assert np.all(np.isclose(resvals, origvals))
-
-        del img
-        img = None
-
-
-def test_image_resample_4d(seed):
-
-    fname = 'test.nii.gz'
-
-    with tempdir():
-
-        make_random_image(fname, (10, 10, 10, 10))
-
-        # resample one volume
-        img = fslimage.Image(fname)
-        slc = (slice(None), slice(None), slice(None), 3)
-        resampled = img.resample(img.shape[:3], slc)[0]
-        assert np.all(resampled == img[..., 3])
-
-        # resample up
-        resampled = img.resample((15, 15, 15), slc)[0]
-        assert tuple(resampled.shape) == (15, 15, 15)
-
-        # resample down
-        resampled = img.resample((5, 5, 5), slc)[0]
-        assert tuple(resampled.shape) == (5, 5, 5)
-
-        # resample the entire image
-        resampled = img.resample((15, 15, 15, 10), None)[0]
-        assert tuple(resampled.shape) == (15, 15, 15, 10)
-
-        resampled = img.resample((5, 5, 5, 10), None)[0]
-        assert tuple(resampled.shape) == (5, 5, 5, 10)
-
-        # resample along the fourth dim
-        resampled = img.resample((15, 15, 15, 15), None)[0]
-        assert tuple(resampled.shape) == (15, 15, 15, 15)
-
-        resampled = img.resample((5, 5, 5, 15), None)[0]
-        assert tuple(resampled.shape) == (5, 5, 5, 15)
-
-        del img
-        del resampled
-        img = None
-        resampled = None
-
-
 def  test_Image_init_xform_nifti1():  _test_Image_init_xform(1)
 def  test_Image_init_xform_nifti2():  _test_Image_init_xform(2)
 def _test_Image_init_xform(imgtype):
diff --git a/tests/test_image_resample.py b/tests/test_image_resample.py
new file mode 100644
index 0000000000000000000000000000000000000000..8802ded1eeee128e0d4cc39f833bb27e213527da
--- /dev/null
+++ b/tests/test_image_resample.py
@@ -0,0 +1,225 @@
+#!/usr/bin/env python
+
+
+import itertools as it
+import numpy     as np
+
+import pytest
+
+import fsl.data.image           as     fslimage
+import fsl.utils.transform      as     transform
+import fsl.utils.image.resample as     resample
+
+from . import make_random_image
+
+
+def test_resample(seed):
+
+    # Random base image shapes
+    for i in range(25):
+
+        shape = np.random.randint(5, 50, 3)
+        img = fslimage.Image(make_random_image(dims=shape))
+
+        # bad shape
+        with pytest.raises(ValueError):
+            resample.resample(img, (10, 10))
+        with pytest.raises(ValueError):
+            resample.resample(img, (10, 10, 10, 10))
+
+        # resampling to the same shape should be a no-op
+        samei, samex = resample.resample(img, shape)
+        assert np.all(samei == img[:])
+        assert np.all(samex == img.voxToWorldMat)
+
+        # Random resampled image shapes
+        for j in range(10):
+
+            rshape        = np.random.randint(5, 50, 3)
+            resampled, xf = resample.resample(img, rshape, order=0)
+
+            assert tuple(resampled.shape) == tuple(rshape)
+
+            # We used nearest neighbour interp, so the
+            # values in the resampled image should match
+            # corresponding values in the original. Let's
+            # check some whynot.
+            restestcoords = np.array([
+                np.random.randint(0, rshape[0], 100),
+                np.random.randint(0, rshape[1], 100),
+                np.random.randint(0, rshape[2], 100)]).T
+
+            resx,  resy,  resz  = restestcoords.T
+            resvals  = resampled[resx, resy, resz]
+
+            res2orig = transform.concat(img.worldToVoxMat, xf)
+
+            origtestcoords = transform.transform(restestcoords, res2orig)
+
+            # remove any coordinates which are out of
+            # bounds in the original image space, or
+            # are right on a voxel boundary (where the
+            # nn interp could have gone either way), or
+            # have value == 0 in the resampled space.
+            out = ((origtestcoords < 0)            |
+                   (origtestcoords >= shape - 0.5) |
+                   (np.isclose(np.modf(origtestcoords)[0], 0.5)))
+            out = np.any(out, axis=1) | (resvals == 0)
+
+            origtestcoords = np.array(origtestcoords.round(), dtype=np.int)
+
+            origtestcoords = origtestcoords[~out, :]
+            restestcoords  = restestcoords[ ~out, :]
+
+            resx,  resy,  resz  = restestcoords.T
+            origx, origy, origz = origtestcoords.T
+
+            origvals = img[:][origx, origy, origz]
+            resvals  = resampled[resx, resy, resz]
+
+            assert np.all(np.isclose(resvals, origvals))
+
+
+def test_resample_4d(seed):
+
+    # resample one volume
+    img = fslimage.Image(make_random_image(dims=(10, 10, 10, 10)))
+    slc = (slice(None), slice(None), slice(None), 3)
+    resampled = resample.resample(img, img.shape[:3], slc)[0]
+    assert np.all(resampled == img[..., 3])
+
+    # resample up
+    resampled = resample.resample(img, (15, 15, 15), slc)[0]
+    assert tuple(resampled.shape) == (15, 15, 15)
+
+    # resample down
+    resampled = resample.resample(img, (5, 5, 5), slc)[0]
+    assert tuple(resampled.shape) == (5, 5, 5)
+
+    # resample the entire image
+    resampled = resample.resample(img, (15, 15, 15, 10), None)[0]
+    assert tuple(resampled.shape) == (15, 15, 15, 10)
+
+    resampled = resample.resample(img, (5, 5, 5, 10), None)[0]
+    assert tuple(resampled.shape) == (5, 5, 5, 10)
+
+    # resample along the fourth dim
+    resampled = resample.resample(img, (15, 15, 15, 15), None)[0]
+    assert tuple(resampled.shape) == (15, 15, 15, 15)
+
+    resampled = resample.resample(img, (5, 5, 5, 15), None)[0]
+    assert tuple(resampled.shape) == (5, 5, 5, 15)
+
+
+def test_resample_origin(seed):
+
+    img = fslimage.Image(make_random_image(dims=(10, 10, 10)))
+
+    # with origin='corner', image
+    # bounding boxes should match
+    for i in range(25):
+        shape = np.random.randint(5, 50, 3)
+        res = resample.resample(img, shape, origin='corner')
+        res = fslimage.Image(res[0], xform=res[1])
+        imgb = transform.axisBounds(img.shape, img.voxToWorldMat)
+        resb = transform.axisBounds(res.shape, res.voxToWorldMat)
+        assert np.all(np.isclose(imgb, resb, rtol=1e-5, atol=1e-5))
+
+    # with origin='centre' image
+    # bounding boxes should be offset
+    # by (size_resampled - size_orig) / 2
+    for i in range(25):
+        shape = np.random.randint(5, 50, 3)
+        res = resample.resample(img, shape, origin='centre')
+        res = fslimage.Image(res[0], xform=res[1])
+        off = (np.array(img.shape) / np.array(res.shape) - 1) / 2
+        imgb = np.array(transform.axisBounds(img.shape, img.voxToWorldMat))
+        resb = np.array(transform.axisBounds(res.shape, res.voxToWorldMat))
+        assert np.all(np.isclose(imgb, resb + off, rtol=1e-5, atol=1e-5))
+
+    # with origin='corner', using
+    # linear interp, when we down-
+    # sample an image to a shape
+    # that divides evenly into the
+    # original shape, a downsampled
+    # voxel should equal the average
+    # of the original voxels inside
+    # it.
+    res = resample.resample(img, (5, 5, 5), smooth=False, origin='corner')[0]
+    for x, y, z in it.product(range(5), range(5), range(5)):
+        block = img[x * 2: x * 2 + 2,
+                    y * 2: y * 2 + 2,
+                    z * 2: z * 2 + 2]
+        assert np.isclose(res[x, y, z], block.mean())
+
+
+def test_resampleToPixdims():
+
+    img          = fslimage.Image(make_random_image(dims=(10, 10, 10)))
+    imglo, imghi = transform.axisBounds(img.shape, img.voxToWorldMat)
+    oldpix       = np.array(img.pixdim, dtype=np.float)
+    oldshape     = np.array(img.shape,  dtype=np.float)
+
+    for i, origin in it.product(range(25), ('centre', 'corner')):
+
+        # random pixdims in the range 0.1 - 5.0
+        newpix   = 0.1 + 4.9 * np.random.random(3)
+        expshape = np.round(oldshape * (oldpix / newpix))
+
+        res = resample.resampleToPixdims(img, newpix, origin=origin)
+        res = fslimage.Image(res[0], xform=res[1])
+        reslo, reshi = transform.axisBounds(res.shape, res.voxToWorldMat)
+        resfov       = reshi - reslo
+        expfov       = newpix * res.shape
+
+        assert np.all(np.isclose(res.shape,  expshape))
+        assert np.all(np.isclose(res.pixdim, newpix))
+        assert np.all(np.isclose(resfov, expfov, rtol=1e-2, atol=1e-2))
+
+        if origin == 'corner':
+            assert np.all(np.isclose(imglo, reslo))
+            assert np.all(np.isclose(reshi, reslo + expfov,
+                                     rtol=1e-2, atol=1e-2))
+
+
+
+def test_resampleToReference():
+
+    def random_v2w():
+        return transform.compose(
+            0.25   + 4.75      * np.random.random(3),
+            -50    + 100       * np.random.random(3),
+            -np.pi + 2 * np.pi * np.random.random(3))
+
+    # Basic test - output has same
+    # dimensions/space as reference
+    for i in range(25):
+
+        ishape = np.random.randint(5, 50, 3)
+        rshape = np.random.randint(5, 50, 3)
+        iv2w   = random_v2w()
+        rv2w   = random_v2w()
+        img    = fslimage.Image(make_random_image(dims=ishape, xform=iv2w))
+        ref    = fslimage.Image(make_random_image(dims=rshape, xform=rv2w))
+        res    = resample.resampleToReference(img, ref)
+        res    = fslimage.Image(res[0], header=ref.header)
+
+        assert res.sameSpace(ref)
+
+    # More specific test - output
+    # data gets transformed correctly
+    # into reference space
+
+    img          = np.zeros((5, 5, 5), dtype=np.float)
+    img[1, 1, 1] = 1
+    img          = fslimage.Image(img)
+
+    refv2w = transform.scaleOffsetXform([1, 1, 1], [-1, -1, -1])
+    ref    = np.zeros((5, 5, 5), dtype=np.float)
+    ref    = fslimage.Image(ref, xform=refv2w)
+    res    = resample.resampleToReference(img, ref, order=0)
+
+    exp          = np.zeros((5, 5, 5), dtype=np.float)
+    exp[2, 2, 2] = 1
+
+    assert np.all(np.isclose(res[0], exp))
diff --git a/tests/test_imglob.py b/tests/test_imglob.py
index 333fbd104c8599dda439543ee9d0d339df126cfd..9a8e838d1321fd45c91c1d94188700067423097f 100644
--- a/tests/test_imglob.py
+++ b/tests/test_imglob.py
@@ -5,14 +5,11 @@
 # Author: Paul McCarthy <pauldmccarthy@gmail.com>
 #
 
-import mock
 import pytest
 
 import fsl.scripts.imglob as imglob
 
-from tests import testdir
-from tests import CaptureStdout
-
+from . import testdir
 
 
 def test_imglob_shouldPass():
@@ -88,60 +85,3 @@ def test_imglob_shouldFail():
 
     with pytest.raises(ValueError):
         imglob.imglob([], 'bag')
-
-
-def test_imglob_script_shouldPass():
-
-    # (files to create, args, expected)
-    tests = [
-        ('file.hdr file.img', 'file',                          'file'),
-        ('file.hdr file.img', '-extension file',               'file.hdr'),
-        ('file.hdr file.img', '-extensions file',              'file.hdr file.img'),
-        ('file.hdr file.img', 'file.hdr',                      'file'),
-        ('file.hdr file.img', ' -extension file.hdr',           'file.hdr'),
-        ('file.hdr file.img', '-extensions file.hdr',          'file.hdr file.img'),
-        ('file.hdr file.img', 'file.img',                      'file'),
-        ('file.hdr file.img', '-extension file.img',           'file.hdr'),
-        ('file.hdr file.img', '-extensions file.img',          'file.hdr file.img'),
-        ('file.hdr file.img', 'file.hdr file.img',             'file'),
-        ('file.hdr file.img', '-extension file.hdr file.img',  'file.hdr'),
-        ('file.hdr file.img', '-extensions file.hdr file.img', 'file.hdr file.img'),
-        ('file.hdr file.img', 'file file.img',                 'file'),
-        ('file.hdr file.img', '-extension file file.img',      'file.hdr'),
-        ('file.hdr file.img', '-extensions file file.img',     'file.hdr file.img'),
-
-        # no file or incomplete prefix
-        ('file.hdr file.img', 'bag',             ''),
-        ('file.hdr file.img', '-extension  bag', ''),
-        ('file.hdr file.img', '-extensions bag', ''),
-        ('file.hdr file.img', 'fi',              ''),
-        ('file.hdr file.img', '-extension  fi',  ''),
-        ('file.hdr file.img', '-extensions fi',  ''),
-    ]
-
-    capture = CaptureStdout()
-
-    for to_create, args, expected in tests:
-        with testdir(to_create.split()) as td:
-
-            capture.reset()
-
-            with capture:
-                assert imglob.main(args.split()) == 0
-
-            assert capture.stdout.strip().split() == expected.split()
-
-
-def test_imglob_script_shouldFail():
-
-    capture = CaptureStdout()
-
-    with capture:
-        assert imglob.main([]) != 0
-
-    assert capture.stdout.strip().lower().startswith('usage:')
-
-    with capture, mock.patch('sys.argv', ['imglob']):
-        assert imglob.main() != 0
-
-    assert capture.stdout.strip().lower().startswith('usage:')
diff --git a/tests/test_immv_imcp.py b/tests/test_immv_imcp.py
index 215b4dc4062245786d4b37af190d626928cefc1c..dd1273d78e162a3ee96638cfba4b12850d1b080f 100644
--- a/tests/test_immv_imcp.py
+++ b/tests/test_immv_imcp.py
@@ -10,26 +10,18 @@ from __future__ import print_function
 
 
 import os.path    as op
-import itertools  as it
-import subprocess as sp
 import               os
 import               shutil
 import               tempfile
 
-import               pytest
-
 import nibabel as nib
 
-from   fsl.utils.tempdir import tempdir
 import fsl.utils.imcp        as imcp
-import fsl.scripts.imcp      as imcp_script
-import fsl.scripts.immv      as immv_script
 import fsl.data.image        as fslimage
 
 from . import make_random_image
 from . import make_dummy_file
 from . import looks_like_image
-from . import cleardir
 
 
 real_print = print
@@ -111,415 +103,6 @@ def checkFilesToExpect(files, outdir, outputType, datahashes):
         checkImageHash(f, h)
 
 
-def test_imcp_script_shouldPass(move=False):
-
-
-    # The imcp/immv scripts should honour the
-    # FSLOUTPUTTYPE env var. If it is unset
-    # or invalid - '' in this case), they
-    # should produce .nii.gz
-    outputTypes = ['NIFTI', 'NIFTI_PAIR', 'NIFTI_GZ', '']
-    reldirs     = ['neutral', 'samedir', 'indir', 'outdir']
-
-    # Test tuples have the following structure (each
-    # item is a string which will be split on spaces):
-    #
-    #   (files_to_create, imcp_args, files_to_expect)
-    #
-    # The files_to_expect is a list of
-    # prefixes - the suffix(es) is(are)
-    # determined by the current outputType
-    #
-    # If multiple files are being copied, the
-    # files_to_create and files_to_expect lists
-    # have to be in the same order.
-    tests = [
-
-        ('a.nii', 'a     b',        'b'),
-        ('a.nii', 'a.nii b',        'b'),
-        ('a.nii', 'a     b.nii',    'b'),
-        ('a.nii', 'a.nii b.nii',    'b'),
-        ('a.nii', 'a     .',        'a'),
-        ('a.nii', 'a.nii .',        'a'),
-        ('a.nii', 'a     b.hdr',    'b'),
-        ('a.nii', 'a     b.img',    'b'),
-        ('a.nii', 'a     b.nii.gz', 'b'),
-        ('a.nii', 'a.nii b.hdr',    'b'),
-        ('a.nii', 'a.nii b.img',    'b'),
-        ('a.nii', 'a.nii b.nii.gz', 'b'),
-
-        ('a.nii.gz', 'a        b',        'b'),
-        ('a.nii.gz', 'a.nii.gz b',        'b'),
-        ('a.nii.gz', 'a        b.nii.gz', 'b'),
-        ('a.nii.gz', 'a.nii.gz b.nii.gz', 'b'),
-        ('a.nii.gz', 'a        .',        'a'),
-        ('a.nii.gz', 'a.nii.gz .',        'a'),
-        ('a.nii.gz', 'a        b.hdr',    'b'),
-        ('a.nii.gz', 'a        b.img',    'b'),
-        ('a.nii.gz', 'a        b.nii',    'b'),
-        ('a.nii.gz', 'a.nii.gz b.hdr',    'b'),
-        ('a.nii.gz', 'a.nii.gz b.img',    'b'),
-        ('a.nii.gz', 'a.nii.gz b.nii',    'b'),
-
-        ('a.img', 'a        b',        'b'),
-        ('a.img', 'a        b.img',    'b'),
-        ('a.img', 'a        b.hdr',    'b'),
-        ('a.img', 'a        .',        'a'),
-        ('a.img', 'a.img    b',        'b'),
-        ('a.img', 'a.img    b.img',    'b'),
-        ('a.img', 'a.img    b.hdr',    'b'),
-        ('a.img', 'a.img    .',        'a'),
-        ('a.img', 'a.hdr    b',        'b'),
-        ('a.img', 'a.hdr    b.img',    'b'),
-        ('a.img', 'a.hdr    b.hdr',    'b'),
-        ('a.img', 'a.hdr    .',        'a'),
-
-        ('a.img', 'a        b.nii',    'b'),
-        ('a.img', 'a        b.nii.gz', 'b'),
-        ('a.img', 'a        .',        'a'),
-        ('a.img', 'a.hdr    b.nii',    'b'),
-        ('a.img', 'a.hdr    b.nii.gz', 'b'),
-        ('a.img', 'a.hdr    .',        'a'),
-        ('a.img', 'a.img    b.nii',    'b'),
-        ('a.img', 'a.img    b.nii.gz', 'b'),
-        ('a.img', 'a.img    .',        'a'),
-
-        ('a.nii b.nii', 'a     b     .',   'a b'),
-        ('a.nii b.nii', 'a     b.nii .',   'a b'),
-        ('a.nii b.nii', 'a.nii b     .',   'a b'),
-        ('a.nii b.nii', 'a.nii b.nii .',   'a b'),
-
-        ('a.img b.img', 'a     b     .',   'a b'),
-        ('a.img b.img', 'a     b.img .',   'a b'),
-        ('a.img b.img', 'a     b.hdr .',   'a b'),
-        ('a.img b.img', 'a.img b     .',   'a b'),
-        ('a.img b.img', 'a.img b.img .',   'a b'),
-        ('a.img b.img', 'a.img b.hdr .',   'a b'),
-        ('a.img b.img', 'a.hdr b     .',   'a b'),
-        ('a.img b.img', 'a.hdr b.img .',   'a b'),
-        ('a.img b.img', 'a.hdr b.hdr .',   'a b'),
-
-        ('a.nii.gz b.nii.gz', 'a        b        .',   'a b'),
-        ('a.nii.gz b.nii.gz', 'a        b.nii.gz .',   'a b'),
-        ('a.nii.gz b.nii.gz', 'a.nii.gz b        .',   'a b'),
-        ('a.nii.gz b.nii.gz', 'a.nii.gz b.nii.gz .',   'a b'),
-
-        # Heterogenous inputs
-        ('a.nii b.nii.gz', 'a     b        .', 'a b'),
-        ('a.nii b.nii.gz', 'a     b.nii.gz .', 'a b'),
-        ('a.nii b.nii.gz', 'a.nii b        .', 'a b'),
-        ('a.nii b.nii.gz', 'a.nii b.nii.gz .', 'a b'),
-        ('a.nii b.img',    'a     b        .', 'a b'),
-        ('a.nii b.img',    'a     b.img    .', 'a b'),
-        ('a.nii b.img',    'a     b.hdr    .', 'a b'),
-        ('a.nii b.img',    'a.nii b        .', 'a b'),
-        ('a.nii b.img',    'a.nii b.img    .', 'a b'),
-        ('a.nii b.img',    'a.nii b.hdr    .', 'a b'),
-        ('a.img b.nii',    'a     b        .', 'a b'),
-        ('a.img b.nii',    'a     b.nii    .', 'a b'),
-        ('a.img b.nii',    'a.img b        .', 'a b'),
-        ('a.img b.nii',    'a.img b.nii    .', 'a b'),
-        ('a.img b.nii',    'a.hdr b        .', 'a b'),
-        ('a.img b.nii',    'a.hdr b.nii    .', 'a b'),
-
-        # Duplicate inputs
-        ('a.img',       'a     a                 .', 'a'),
-        ('a.img',       'a     a.img             .', 'a'),
-        ('a.img',       'a     a.hdr             .', 'a'),
-        ('a.img',       'a.img a                 .', 'a'),
-        ('a.img',       'a.img a.img             .', 'a'),
-        ('a.img',       'a.img a.hdr             .', 'a'),
-        ('a.img',       'a.hdr a                 .', 'a'),
-        ('a.img',       'a.hdr a.img             .', 'a'),
-        ('a.img',       'a.hdr a.hdr             .', 'a'),
-
-        ('a.img b.img', 'a     a     b     b     .', 'a b'),
-        ('a.img b.img', 'a     a     b     b.img .', 'a b'),
-        ('a.img b.img', 'a     a     b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a     a     b.img b     .', 'a b'),
-        ('a.img b.img', 'a     a     b.img b.img .', 'a b'),
-        ('a.img b.img', 'a     a     b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a     a.img b     b     .', 'a b'),
-        ('a.img b.img', 'a     a.img b     b.img .', 'a b'),
-        ('a.img b.img', 'a     a.img b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a     a.img b.img b     .', 'a b'),
-        ('a.img b.img', 'a     a.img b.img b.img .', 'a b'),
-        ('a.img b.img', 'a     a.img b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b     b     .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b     b.img .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b.img b     .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b.img b.img .', 'a b'),
-        ('a.img b.img', 'a     a.hdr b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a     b     b     .', 'a b'),
-        ('a.img b.img', 'a.img a     b     b.img .', 'a b'),
-        ('a.img b.img', 'a.img a     b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a     b.img b     .', 'a b'),
-        ('a.img b.img', 'a.img a     b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.img a     b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a.img b     b     .', 'a b'),
-        ('a.img b.img', 'a.img a.img b     b.img .', 'a b'),
-        ('a.img b.img', 'a.img a.img b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a.img b.img b     .', 'a b'),
-        ('a.img b.img', 'a.img a.img b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.img a.img b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b     b     .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b     b.img .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b.img b     .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.img a.hdr b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b     b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b     b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b.img b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a     b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b     b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b     b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b.img b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a.img b.img b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b     b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b     b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b     b.hdr .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b.img b     .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b.img b.img .', 'a b'),
-        ('a.img b.img', 'a.hdr a.hdr b.img b.hdr .', 'a b'),
-
-        # Inputs which cause the destination
-        # to be overwritten - this should be
-        # ok, the destination should be the
-        # last specified input. The order of
-        # files_to_create has to match the
-        # imcp_args order, otherwise my bodgy
-        # checkFilesToExpect function will
-        # break.
-        ('a.nii    a.img',             'a.nii    a.img             .', 'a'),
-        ('a.img    a.nii',             'a.img    a.nii             .', 'a'),
-        ('a.nii    a.img    a.nii.gz', 'a.nii    a.img    a.nii.gz .', 'a'),
-        ('a.nii    a.nii.gz a.img   ', 'a.nii    a.nii.gz a.img    .', 'a'),
-        ('a.img    a.nii.gz a.nii   ', 'a.img    a.nii.gz a.nii    .', 'a'),
-        ('a.img    a.nii    a.nii.gz', 'a.img    a.nii    a.nii.gz .', 'a'),
-        ('a.nii.gz a.img    a.nii   ', 'a.nii.gz a.img    a.nii    .', 'a'),
-        ('a.nii.gz a.nii    a.img   ', 'a.nii.gz a.nii    a.img    .', 'a'),
-    ]
-
-    indir    = tempfile.mkdtemp()
-    outdir   = tempfile.mkdtemp()
-    startdir = os.getcwd()
-
-    try:
-
-        for outputType, reldir in it.product(outputTypes, reldirs):
-
-            os.environ['FSLOUTPUTTYPE'] = outputType
-
-            for files_to_create, imcp_args, files_to_expect in tests:
-
-                imageHashes = []
-
-                print()
-                print('files_to_create: ', files_to_create)
-                print('imcp_args:       ', imcp_args)
-                print('files_to_expect: ', files_to_expect)
-
-                for i, fname in enumerate(files_to_create.split()):
-                    imageHashes.append(makeImage(op.join(indir, fname)))
-
-                imcp_args = imcp_args.split()
-
-                tindir  = indir
-                toutdir = outdir
-
-                if   reldir == 'neutral': reldir = startdir
-                elif reldir == 'indir':   reldir = tindir
-                elif reldir == 'outdir':  reldir = toutdir
-                elif reldir == 'samedir':
-                    reldir  = tindir
-                    toutdir = tindir
-
-                    if not move:
-
-                        infiles = os.listdir(tindir)
-
-                        files_to_expect = files_to_expect +  ' ' + \
-                                          ' '.join(infiles)
-
-                        for inf in infiles:
-                            img     = nib.load(op.join(tindir, inf),
-                                               mmap=False)
-                            imghash = hash(img.get_data().tobytes())
-                            img = None
-                            imageHashes.append(imghash)
-
-                print('adj files_to_expect: ', files_to_expect)
-
-                os.chdir(reldir)
-
-                imcp_args[:-1] = [op.join(tindir, a) for a in imcp_args[:-1]]
-                imcp_args[ -1] =  op.join(toutdir, imcp_args[-1])
-
-                for i, a in enumerate(imcp_args):
-                    if op.splitdrive(a)[0] == op.splitdrive(reldir)[0]:
-                        imcp_args[i] = op.relpath(a, reldir)
-
-                print('indir before:    ', os.listdir(tindir))
-                print('outdir before:   ', os.listdir(toutdir))
-
-                if move: result = immv_script.main(imcp_args)
-                else:    result = imcp_script.main(imcp_args)
-
-                print('indir after:     ', os.listdir(tindir))
-                print('outdir after:    ', os.listdir(toutdir))
-
-                assert result == 0
-
-                checkFilesToExpect(
-                    files_to_expect, toutdir, outputType, imageHashes)
-
-                # too hard if indir == outdir
-                if move and tindir != toutdir:
-                    infiles = os.listdir(tindir)
-                    infiles = [f for f in infiles if op.isfile(f)]
-                    infiles = [f for f in infiles if op.isfile(f)]
-                    assert len(infiles) == 0
-
-                cleardir(indir)
-                cleardir(outdir)
-
-    finally:
-        os.chdir(startdir)
-        shutil.rmtree(indir)
-        shutil.rmtree(outdir)
-
-
-@pytest.mark.noroottest
-def test_imcp_script_shouldFail(move=False):
-
-    # - len(srcs) > 1 and dest is not dir
-    # - input not readable
-    # - move=True and input not deleteable
-    # - ambiguous inputs
-    # - input is incomplete pair (e.g. .hdr
-    #   without corresponding .img)
-
-    # a.img
-    # a.hdr
-    # a.nii
-    #
-    # FAIL: imcp a           dest
-
-    # (files_to_create, imcp_args[, preproc])
-    tests = [
-
-        # non-existent input
-        ('',      'a     b'),
-        ('a.img', 'b.img a'),
-
-        # dest is non-existent dir
-        ('a.img',       'a.img non_existent_dir/'),
-
-        # len(srcs) > 1, but dest is not dir
-        ('a.img b.img', 'a b c.img'),
-
-        # destination not writeable
-        ('a.img b.img', 'a b ./readonly_dir', 'mkdir outdir/readonly_dir; chmod a-wx outdir/readonly_dir'),
-        ('a.img',       'a    b',             'mkdir outdir/b.nii.gz; chmod a-wx outdir/b.nii.gz'),
-
-        # input not readable
-        ('a.img', 'a b', 'chmod a-rwx indir/a.img'),
-
-        # ambiguous input
-        ('a.img a.nii', 'a b'),
-
-        # input is part of incomplete pair
-        ('a.img', 'a     b', 'rm indir/a.hdr'),
-        ('a.img', 'a.img b', 'rm indir/a.hdr'),
-        ('a.img', 'a.hdr b', 'rm indir/a.hdr'),
-        ('a.img', 'a     b', 'rm indir/a.img'),
-        ('a.img', 'a.img b', 'rm indir/a.img'),
-        ('a.img', 'a.hdr b', 'rm indir/a.img'),
-    ]
-
-    if move:
-        tests = tests + [
-            # Input not deletable
-            ('a.img', 'a b', 'chmod a-wx indir'),
-            ('a.img', 'a b', 'chmod a-wx indir'),
-            ('a.nii', 'a b', 'chmod a-wx indir')
-        ]
-
-    indir  = tempfile.mkdtemp()
-    outdir = tempfile.mkdtemp()
-
-    try:
-        for test in tests:
-
-            files_to_create = test[0]
-            imcp_args       = test[1]
-
-            if len(test) == 3: preproc = test[2]
-            else:              preproc = None
-
-            files_to_create = files_to_create.split()
-            imcp_args       = imcp_args      .split()
-
-            for fname in files_to_create:
-                makeImage(op.join(indir, fname))
-
-            imcp_args[:-1] = [op.join(indir, a) for a in imcp_args[:-1]]
-            imcp_args[ -1] =  op.join(outdir, imcp_args[-1])
-
-            if preproc is not None:
-                for cmd in preproc.split(';'):
-                    cmd = cmd.replace('indir', indir).replace('outdir', outdir)
-                    sp.call(cmd.split())
-
-            print('calling {} {}'.format('immv' if move else 'imcp',
-                                         ' '.join(imcp_args)))
-
-            print('indir before:   {}'.format(os.listdir(indir)))
-            print('out dir before: {}'.format(os.listdir(outdir)))
-
-            if move: result = immv_script.main(imcp_args)
-            else:    result = imcp_script.main(imcp_args)
-
-            print('indir after:   {}'.format(os.listdir(indir)))
-            print('out dir after: {}'.format(os.listdir(outdir)))
-
-            assert result != 0
-
-            sp.call('chmod u+rwx {}'.format(indir) .split())
-            sp.call('chmod u+rwx {}'.format(outdir).split())
-
-            cleardir(indir)
-            cleardir(outdir)
-
-    finally:
-        shutil.rmtree(indir)
-        shutil.rmtree(outdir)
-
-    # Other failure cases
-    if move: assert immv_script.main()       != 0
-    else:    assert imcp_script.main()       != 0
-    if move: assert immv_script.main([])     != 0
-    else:    assert imcp_script.main([])     != 0
-    if move: assert immv_script.main(['wa']) != 0
-    else:    assert imcp_script.main(['wa']) != 0
-
-
-
-def test_immv_script_shouldPass():
-    test_imcp_script_shouldPass(move=True)
-
-
-@pytest.mark.noroottest
-def test_immv_script_shouldFail():
-    test_imcp_script_shouldFail(move=True)
-
-
 def test_imcp_shouldPass(move=False):
 
 
@@ -731,27 +314,3 @@ def test_imcp_shouldPass(move=False):
 
 def test_immv_shouldPass():
     test_imcp_shouldPass(move=True)
-
-def test_imcp_badExt():
-    with tempdir():
-
-        with open('file.nii.gz', 'wt') as f:
-            f.write('1')
-
-        result = imcp_script.main(['file.nii', 'dest'])
-
-        assert result == 0
-        assert op.exists('dest.nii.gz')
-
-
-
-def test_immv_badExt():
-    with tempdir():
-
-        with open('file.nii.gz', 'wt') as f:
-            f.write('1')
-
-        result = immv_script.main(['file.nii', 'dest'])
-
-        assert result == 0
-        assert op.exists('dest.nii.gz')
diff --git a/tests/test_scripts/__init__.py b/tests/test_scripts/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/tests/test_atlasq_list_summary.py b/tests/test_scripts/test_atlasq_list_summary.py
similarity index 98%
rename from tests/test_atlasq_list_summary.py
rename to tests/test_scripts/test_atlasq_list_summary.py
index 310ef9212af987dd96f81b6340a9c8307d7031fa..6fbebb7b6ee19d9feb31e71e98923d18c03a2c08 100644
--- a/tests/test_atlasq_list_summary.py
+++ b/tests/test_scripts/test_atlasq_list_summary.py
@@ -13,7 +13,7 @@ import              pytest
 import fsl.data.atlases   as fslatlases
 import fsl.scripts.atlasq as fslatlasq
 
-from . import CaptureStdout
+from .. import CaptureStdout
 
 
 pytestmark = pytest.mark.fsltest
diff --git a/tests/test_atlasq_ohi.py b/tests/test_scripts/test_atlasq_ohi.py
similarity index 98%
rename from tests/test_atlasq_ohi.py
rename to tests/test_scripts/test_atlasq_ohi.py
index 2c0e4e7aadbeea076278f215b56a5e851458f43d..c2f2348a532c04f9819c28d980393bc3857c7746 100644
--- a/tests/test_atlasq_ohi.py
+++ b/tests/test_scripts/test_atlasq_ohi.py
@@ -18,9 +18,9 @@ import numpy as np
 import fsl.scripts.atlasq as fslatlasq
 import fsl.data.atlases   as fslatlases
 
-from . import (tempdir,
-               make_random_mask,
-               CaptureStdout)
+from .. import (tempdir,
+                make_random_mask,
+                CaptureStdout)
 
 
 pytestmark = pytest.mark.fsltest
diff --git a/tests/test_atlasq_query.py b/tests/test_scripts/test_atlasq_query.py
similarity index 98%
rename from tests/test_atlasq_query.py
rename to tests/test_scripts/test_atlasq_query.py
index 987cd01141211a28b26739bdb524d8bee3d9f039..a36e1eb5fa4f086b0e28bcd4440f08958ebb3d05 100644
--- a/tests/test_atlasq_query.py
+++ b/tests/test_scripts/test_atlasq_query.py
@@ -16,13 +16,14 @@ import scipy.ndimage as ndi
 import                  pytest
 
 import fsl.utils.transform as transform
+import fsl.utils.image.resample as resample
 import fsl.data.atlases    as fslatlases
 import fsl.data.image      as fslimage
 import fsl.scripts.atlasq  as fslatlasq
 
-from . import (tempdir,
-               make_random_mask,
-               CaptureStdout)
+from .. import (tempdir,
+                make_random_mask,
+                CaptureStdout)
 
 
 pytestmark = pytest.mark.fsltest
@@ -329,7 +330,8 @@ def _gen_mask_query(atlas, use_label, q_type, q_in, res):
             # resampled into the atlas resolution,
             # it is still either in or out of the
             # atlas space
-            mask, xform = mask.resample(a.shape[:3], dtype=np.float32, order=1)
+            mask, xform = resample.resample(
+                mask, a.shape[:3], dtype=np.float32, order=1)
 
             thres = np.percentile(mask[mask > 0], 75)
 
diff --git a/tests/test_fsl_ents.py b/tests/test_scripts/test_fsl_ents.py
similarity index 100%
rename from tests/test_fsl_ents.py
rename to tests/test_scripts/test_fsl_ents.py
diff --git a/tests/test_scripts/test_imglob.py b/tests/test_scripts/test_imglob.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f8e934509ce68359d17a8a8794bbaaa47f1a691
--- /dev/null
+++ b/tests/test_scripts/test_imglob.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+
+import mock
+
+import fsl.scripts.imglob as imglob
+
+from .. import testdir
+from .. import CaptureStdout
+
+
+def test_imglob_script_shouldPass():
+
+    # (files to create, args, expected)
+    tests = [
+        ('file.hdr file.img', 'file',                          'file'),
+        ('file.hdr file.img', '-extension file',               'file.hdr'),
+        ('file.hdr file.img', '-extensions file',              'file.hdr file.img'),
+        ('file.hdr file.img', 'file.hdr',                      'file'),
+        ('file.hdr file.img', ' -extension file.hdr',           'file.hdr'),
+        ('file.hdr file.img', '-extensions file.hdr',          'file.hdr file.img'),
+        ('file.hdr file.img', 'file.img',                      'file'),
+        ('file.hdr file.img', '-extension file.img',           'file.hdr'),
+        ('file.hdr file.img', '-extensions file.img',          'file.hdr file.img'),
+        ('file.hdr file.img', 'file.hdr file.img',             'file'),
+        ('file.hdr file.img', '-extension file.hdr file.img',  'file.hdr'),
+        ('file.hdr file.img', '-extensions file.hdr file.img', 'file.hdr file.img'),
+        ('file.hdr file.img', 'file file.img',                 'file'),
+        ('file.hdr file.img', '-extension file file.img',      'file.hdr'),
+        ('file.hdr file.img', '-extensions file file.img',     'file.hdr file.img'),
+
+        # no file or incomplete prefix
+        ('file.hdr file.img', 'bag',             ''),
+        ('file.hdr file.img', '-extension  bag', ''),
+        ('file.hdr file.img', '-extensions bag', ''),
+        ('file.hdr file.img', 'fi',              ''),
+        ('file.hdr file.img', '-extension  fi',  ''),
+        ('file.hdr file.img', '-extensions fi',  ''),
+    ]
+
+    capture = CaptureStdout()
+
+    for to_create, args, expected in tests:
+        with testdir(to_create.split()) as td:
+
+            capture.reset()
+
+            with capture:
+                assert imglob.main(args.split()) == 0
+
+            assert capture.stdout.strip().split() == expected.split()
+
+
+def test_imglob_script_shouldFail():
+
+    capture = CaptureStdout()
+
+    with capture:
+        assert imglob.main([]) != 0
+
+    assert capture.stdout.strip().lower().startswith('usage:')
+
+    with capture, mock.patch('sys.argv', ['imglob']):
+        assert imglob.main() != 0
+
+    assert capture.stdout.strip().lower().startswith('usage:')
diff --git a/tests/test_scripts/test_immv_imcp.py b/tests/test_scripts/test_immv_imcp.py
new file mode 100644
index 0000000000000000000000000000000000000000..607a7da5d761d69c824f5d85f27751b8bf5a395d
--- /dev/null
+++ b/tests/test_scripts/test_immv_imcp.py
@@ -0,0 +1,465 @@
+#!/usr/bin/env python
+
+
+
+import os.path    as op
+import itertools  as it
+import subprocess as sp
+import               os
+import               shutil
+import               tempfile
+
+import               pytest
+
+import nibabel as nib
+
+from   fsl.utils.tempdir import tempdir
+import fsl.scripts.imcp      as imcp_script
+import fsl.scripts.immv      as immv_script
+
+from .. import cleardir
+from ..test_immv_imcp import makeImage, checkImageHash, checkFilesToExpect
+
+
+real_print = print
+
+def print(*args, **kwargs):
+    pass
+
+
+
+
+
+def test_imcp_script_shouldPass(move=False):
+
+
+    # The imcp/immv scripts should honour the
+    # FSLOUTPUTTYPE env var. If it is unset
+    # or invalid - '' in this case), they
+    # should produce .nii.gz
+    outputTypes = ['NIFTI', 'NIFTI_PAIR', 'NIFTI_GZ', '']
+    reldirs     = ['neutral', 'samedir', 'indir', 'outdir']
+
+    # Test tuples have the following structure (each
+    # item is a string which will be split on spaces):
+    #
+    #   (files_to_create, imcp_args, files_to_expect)
+    #
+    # The files_to_expect is a list of
+    # prefixes - the suffix(es) is(are)
+    # determined by the current outputType
+    #
+    # If multiple files are being copied, the
+    # files_to_create and files_to_expect lists
+    # have to be in the same order.
+    tests = [
+
+        ('a.nii', 'a     b',        'b'),
+        ('a.nii', 'a.nii b',        'b'),
+        ('a.nii', 'a     b.nii',    'b'),
+        ('a.nii', 'a.nii b.nii',    'b'),
+        ('a.nii', 'a     .',        'a'),
+        ('a.nii', 'a.nii .',        'a'),
+        ('a.nii', 'a     b.hdr',    'b'),
+        ('a.nii', 'a     b.img',    'b'),
+        ('a.nii', 'a     b.nii.gz', 'b'),
+        ('a.nii', 'a.nii b.hdr',    'b'),
+        ('a.nii', 'a.nii b.img',    'b'),
+        ('a.nii', 'a.nii b.nii.gz', 'b'),
+
+        ('a.nii.gz', 'a        b',        'b'),
+        ('a.nii.gz', 'a.nii.gz b',        'b'),
+        ('a.nii.gz', 'a        b.nii.gz', 'b'),
+        ('a.nii.gz', 'a.nii.gz b.nii.gz', 'b'),
+        ('a.nii.gz', 'a        .',        'a'),
+        ('a.nii.gz', 'a.nii.gz .',        'a'),
+        ('a.nii.gz', 'a        b.hdr',    'b'),
+        ('a.nii.gz', 'a        b.img',    'b'),
+        ('a.nii.gz', 'a        b.nii',    'b'),
+        ('a.nii.gz', 'a.nii.gz b.hdr',    'b'),
+        ('a.nii.gz', 'a.nii.gz b.img',    'b'),
+        ('a.nii.gz', 'a.nii.gz b.nii',    'b'),
+
+        ('a.img', 'a        b',        'b'),
+        ('a.img', 'a        b.img',    'b'),
+        ('a.img', 'a        b.hdr',    'b'),
+        ('a.img', 'a        .',        'a'),
+        ('a.img', 'a.img    b',        'b'),
+        ('a.img', 'a.img    b.img',    'b'),
+        ('a.img', 'a.img    b.hdr',    'b'),
+        ('a.img', 'a.img    .',        'a'),
+        ('a.img', 'a.hdr    b',        'b'),
+        ('a.img', 'a.hdr    b.img',    'b'),
+        ('a.img', 'a.hdr    b.hdr',    'b'),
+        ('a.img', 'a.hdr    .',        'a'),
+
+        ('a.img', 'a        b.nii',    'b'),
+        ('a.img', 'a        b.nii.gz', 'b'),
+        ('a.img', 'a        .',        'a'),
+        ('a.img', 'a.hdr    b.nii',    'b'),
+        ('a.img', 'a.hdr    b.nii.gz', 'b'),
+        ('a.img', 'a.hdr    .',        'a'),
+        ('a.img', 'a.img    b.nii',    'b'),
+        ('a.img', 'a.img    b.nii.gz', 'b'),
+        ('a.img', 'a.img    .',        'a'),
+
+        ('a.nii b.nii', 'a     b     .',   'a b'),
+        ('a.nii b.nii', 'a     b.nii .',   'a b'),
+        ('a.nii b.nii', 'a.nii b     .',   'a b'),
+        ('a.nii b.nii', 'a.nii b.nii .',   'a b'),
+
+        ('a.img b.img', 'a     b     .',   'a b'),
+        ('a.img b.img', 'a     b.img .',   'a b'),
+        ('a.img b.img', 'a     b.hdr .',   'a b'),
+        ('a.img b.img', 'a.img b     .',   'a b'),
+        ('a.img b.img', 'a.img b.img .',   'a b'),
+        ('a.img b.img', 'a.img b.hdr .',   'a b'),
+        ('a.img b.img', 'a.hdr b     .',   'a b'),
+        ('a.img b.img', 'a.hdr b.img .',   'a b'),
+        ('a.img b.img', 'a.hdr b.hdr .',   'a b'),
+
+        ('a.nii.gz b.nii.gz', 'a        b        .',   'a b'),
+        ('a.nii.gz b.nii.gz', 'a        b.nii.gz .',   'a b'),
+        ('a.nii.gz b.nii.gz', 'a.nii.gz b        .',   'a b'),
+        ('a.nii.gz b.nii.gz', 'a.nii.gz b.nii.gz .',   'a b'),
+
+        # Heterogenous inputs
+        ('a.nii b.nii.gz', 'a     b        .', 'a b'),
+        ('a.nii b.nii.gz', 'a     b.nii.gz .', 'a b'),
+        ('a.nii b.nii.gz', 'a.nii b        .', 'a b'),
+        ('a.nii b.nii.gz', 'a.nii b.nii.gz .', 'a b'),
+        ('a.nii b.img',    'a     b        .', 'a b'),
+        ('a.nii b.img',    'a     b.img    .', 'a b'),
+        ('a.nii b.img',    'a     b.hdr    .', 'a b'),
+        ('a.nii b.img',    'a.nii b        .', 'a b'),
+        ('a.nii b.img',    'a.nii b.img    .', 'a b'),
+        ('a.nii b.img',    'a.nii b.hdr    .', 'a b'),
+        ('a.img b.nii',    'a     b        .', 'a b'),
+        ('a.img b.nii',    'a     b.nii    .', 'a b'),
+        ('a.img b.nii',    'a.img b        .', 'a b'),
+        ('a.img b.nii',    'a.img b.nii    .', 'a b'),
+        ('a.img b.nii',    'a.hdr b        .', 'a b'),
+        ('a.img b.nii',    'a.hdr b.nii    .', 'a b'),
+
+        # Duplicate inputs
+        ('a.img',       'a     a                 .', 'a'),
+        ('a.img',       'a     a.img             .', 'a'),
+        ('a.img',       'a     a.hdr             .', 'a'),
+        ('a.img',       'a.img a                 .', 'a'),
+        ('a.img',       'a.img a.img             .', 'a'),
+        ('a.img',       'a.img a.hdr             .', 'a'),
+        ('a.img',       'a.hdr a                 .', 'a'),
+        ('a.img',       'a.hdr a.img             .', 'a'),
+        ('a.img',       'a.hdr a.hdr             .', 'a'),
+
+        ('a.img b.img', 'a     a     b     b     .', 'a b'),
+        ('a.img b.img', 'a     a     b     b.img .', 'a b'),
+        ('a.img b.img', 'a     a     b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a     a     b.img b     .', 'a b'),
+        ('a.img b.img', 'a     a     b.img b.img .', 'a b'),
+        ('a.img b.img', 'a     a     b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a     a.img b     b     .', 'a b'),
+        ('a.img b.img', 'a     a.img b     b.img .', 'a b'),
+        ('a.img b.img', 'a     a.img b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a     a.img b.img b     .', 'a b'),
+        ('a.img b.img', 'a     a.img b.img b.img .', 'a b'),
+        ('a.img b.img', 'a     a.img b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b     b     .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b     b.img .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b.img b     .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b.img b.img .', 'a b'),
+        ('a.img b.img', 'a     a.hdr b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a     b     b     .', 'a b'),
+        ('a.img b.img', 'a.img a     b     b.img .', 'a b'),
+        ('a.img b.img', 'a.img a     b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a     b.img b     .', 'a b'),
+        ('a.img b.img', 'a.img a     b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.img a     b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a.img b     b     .', 'a b'),
+        ('a.img b.img', 'a.img a.img b     b.img .', 'a b'),
+        ('a.img b.img', 'a.img a.img b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a.img b.img b     .', 'a b'),
+        ('a.img b.img', 'a.img a.img b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.img a.img b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b     b     .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b     b.img .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b.img b     .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.img a.hdr b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b     b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b     b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b.img b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a     b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b     b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b     b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b.img b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a.img b.img b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b     b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b     b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b     b.hdr .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b.img b     .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b.img b.img .', 'a b'),
+        ('a.img b.img', 'a.hdr a.hdr b.img b.hdr .', 'a b'),
+
+        # Inputs which cause the destination
+        # to be overwritten - this should be
+        # ok, the destination should be the
+        # last specified input. The order of
+        # files_to_create has to match the
+        # imcp_args order, otherwise my bodgy
+        # checkFilesToExpect function will
+        # break.
+        ('a.nii    a.img',             'a.nii    a.img             .', 'a'),
+        ('a.img    a.nii',             'a.img    a.nii             .', 'a'),
+        ('a.nii    a.img    a.nii.gz', 'a.nii    a.img    a.nii.gz .', 'a'),
+        ('a.nii    a.nii.gz a.img   ', 'a.nii    a.nii.gz a.img    .', 'a'),
+        ('a.img    a.nii.gz a.nii   ', 'a.img    a.nii.gz a.nii    .', 'a'),
+        ('a.img    a.nii    a.nii.gz', 'a.img    a.nii    a.nii.gz .', 'a'),
+        ('a.nii.gz a.img    a.nii   ', 'a.nii.gz a.img    a.nii    .', 'a'),
+        ('a.nii.gz a.nii    a.img   ', 'a.nii.gz a.nii    a.img    .', 'a'),
+    ]
+
+    indir    = tempfile.mkdtemp()
+    outdir   = tempfile.mkdtemp()
+    startdir = os.getcwd()
+
+    try:
+
+        for outputType, reldir in it.product(outputTypes, reldirs):
+
+            os.environ['FSLOUTPUTTYPE'] = outputType
+
+            for files_to_create, imcp_args, files_to_expect in tests:
+
+                imageHashes = []
+
+                print()
+                print('files_to_create: ', files_to_create)
+                print('imcp_args:       ', imcp_args)
+                print('files_to_expect: ', files_to_expect)
+
+                for i, fname in enumerate(files_to_create.split()):
+                    imageHashes.append(makeImage(op.join(indir, fname)))
+
+                imcp_args = imcp_args.split()
+
+                tindir  = indir
+                toutdir = outdir
+
+                if   reldir == 'neutral': reldir = startdir
+                elif reldir == 'indir':   reldir = tindir
+                elif reldir == 'outdir':  reldir = toutdir
+                elif reldir == 'samedir':
+                    reldir  = tindir
+                    toutdir = tindir
+
+                    if not move:
+
+                        infiles = os.listdir(tindir)
+
+                        files_to_expect = files_to_expect +  ' ' + \
+                                          ' '.join(infiles)
+
+                        for inf in infiles:
+                            img     = nib.load(op.join(tindir, inf),
+                                               mmap=False)
+                            imghash = hash(img.get_data().tobytes())
+                            img = None
+                            imageHashes.append(imghash)
+
+                print('adj files_to_expect: ', files_to_expect)
+
+                os.chdir(reldir)
+
+                imcp_args[:-1] = [op.join(tindir, a) for a in imcp_args[:-1]]
+                imcp_args[ -1] =  op.join(toutdir, imcp_args[-1])
+
+                for i, a in enumerate(imcp_args):
+                    if op.splitdrive(a)[0] == op.splitdrive(reldir)[0]:
+                        imcp_args[i] = op.relpath(a, reldir)
+
+                print('indir before:    ', os.listdir(tindir))
+                print('outdir before:   ', os.listdir(toutdir))
+
+                if move: result = immv_script.main(imcp_args)
+                else:    result = imcp_script.main(imcp_args)
+
+                print('indir after:     ', os.listdir(tindir))
+                print('outdir after:    ', os.listdir(toutdir))
+
+                assert result == 0
+
+                checkFilesToExpect(
+                    files_to_expect, toutdir, outputType, imageHashes)
+
+                # too hard if indir == outdir
+                if move and tindir != toutdir:
+                    infiles = os.listdir(tindir)
+                    infiles = [f for f in infiles if op.isfile(f)]
+                    infiles = [f for f in infiles if op.isfile(f)]
+                    assert len(infiles) == 0
+
+                cleardir(indir)
+                cleardir(outdir)
+
+    finally:
+        os.chdir(startdir)
+        shutil.rmtree(indir)
+        shutil.rmtree(outdir)
+
+
+@pytest.mark.noroottest
+def test_imcp_script_shouldFail(move=False):
+
+    # - len(srcs) > 1 and dest is not dir
+    # - input not readable
+    # - move=True and input not deleteable
+    # - ambiguous inputs
+    # - input is incomplete pair (e.g. .hdr
+    #   without corresponding .img)
+
+    # a.img
+    # a.hdr
+    # a.nii
+    #
+    # FAIL: imcp a           dest
+
+    # (files_to_create, imcp_args[, preproc])
+    tests = [
+
+        # non-existent input
+        ('',      'a     b'),
+        ('a.img', 'b.img a'),
+
+        # dest is non-existent dir
+        ('a.img',       'a.img non_existent_dir/'),
+
+        # len(srcs) > 1, but dest is not dir
+        ('a.img b.img', 'a b c.img'),
+
+        # destination not writeable
+        ('a.img b.img', 'a b ./readonly_dir', 'mkdir outdir/readonly_dir; chmod a-wx outdir/readonly_dir'),
+        ('a.img',       'a    b',             'mkdir outdir/b.nii.gz; chmod a-wx outdir/b.nii.gz'),
+
+        # input not readable
+        ('a.img', 'a b', 'chmod a-rwx indir/a.img'),
+
+        # ambiguous input
+        ('a.img a.nii', 'a b'),
+
+        # input is part of incomplete pair
+        ('a.img', 'a     b', 'rm indir/a.hdr'),
+        ('a.img', 'a.img b', 'rm indir/a.hdr'),
+        ('a.img', 'a.hdr b', 'rm indir/a.hdr'),
+        ('a.img', 'a     b', 'rm indir/a.img'),
+        ('a.img', 'a.img b', 'rm indir/a.img'),
+        ('a.img', 'a.hdr b', 'rm indir/a.img'),
+    ]
+
+    if move:
+        tests = tests + [
+            # Input not deletable
+            ('a.img', 'a b', 'chmod a-wx indir'),
+            ('a.img', 'a b', 'chmod a-wx indir'),
+            ('a.nii', 'a b', 'chmod a-wx indir')
+        ]
+
+    indir  = tempfile.mkdtemp()
+    outdir = tempfile.mkdtemp()
+
+    try:
+        for test in tests:
+
+            files_to_create = test[0]
+            imcp_args       = test[1]
+
+            if len(test) == 3: preproc = test[2]
+            else:              preproc = None
+
+            files_to_create = files_to_create.split()
+            imcp_args       = imcp_args      .split()
+
+            for fname in files_to_create:
+                makeImage(op.join(indir, fname))
+
+            imcp_args[:-1] = [op.join(indir, a) for a in imcp_args[:-1]]
+            imcp_args[ -1] =  op.join(outdir, imcp_args[-1])
+
+            if preproc is not None:
+                for cmd in preproc.split(';'):
+                    cmd = cmd.replace('indir', indir).replace('outdir', outdir)
+                    sp.call(cmd.split())
+
+            print('calling {} {}'.format('immv' if move else 'imcp',
+                                         ' '.join(imcp_args)))
+
+            print('indir before:   {}'.format(os.listdir(indir)))
+            print('out dir before: {}'.format(os.listdir(outdir)))
+
+            if move: result = immv_script.main(imcp_args)
+            else:    result = imcp_script.main(imcp_args)
+
+            print('indir after:   {}'.format(os.listdir(indir)))
+            print('out dir after: {}'.format(os.listdir(outdir)))
+
+            assert result != 0
+
+            sp.call('chmod u+rwx {}'.format(indir) .split())
+            sp.call('chmod u+rwx {}'.format(outdir).split())
+
+            cleardir(indir)
+            cleardir(outdir)
+
+    finally:
+        shutil.rmtree(indir)
+        shutil.rmtree(outdir)
+
+    # Other failure cases
+    if move: assert immv_script.main()       != 0
+    else:    assert imcp_script.main()       != 0
+    if move: assert immv_script.main([])     != 0
+    else:    assert imcp_script.main([])     != 0
+    if move: assert immv_script.main(['wa']) != 0
+    else:    assert imcp_script.main(['wa']) != 0
+
+
+
+def test_immv_script_shouldPass():
+    test_imcp_script_shouldPass(move=True)
+
+
+@pytest.mark.noroottest
+def test_immv_script_shouldFail():
+    test_imcp_script_shouldFail(move=True)
+
+
+
+def test_imcp_badExt():
+    with tempdir():
+
+        with open('file.nii.gz', 'wt') as f:
+            f.write('1')
+
+        result = imcp_script.main(['file.nii', 'dest'])
+
+        assert result == 0
+        assert op.exists('dest.nii.gz')
+
+
+
+def test_immv_badExt():
+    with tempdir():
+
+        with open('file.nii.gz', 'wt') as f:
+            f.write('1')
+
+        result = immv_script.main(['file.nii', 'dest'])
+
+        assert result == 0
+        assert op.exists('dest.nii.gz')
diff --git a/tests/test_scripts/test_resample_image.py b/tests/test_scripts/test_resample_image.py
new file mode 100644
index 0000000000000000000000000000000000000000..7fc0cf9765fa4dbba615f88e20e5068b6786e8fd
--- /dev/null
+++ b/tests/test_scripts/test_resample_image.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+
+
+import numpy as np
+
+import pytest
+
+import fsl.scripts.resample_image as resample_image
+
+
+import  fsl.utils.transform as transform
+from fsl.utils.tempdir import tempdir
+from fsl.data.image    import Image
+
+from .. import make_random_image
+
+
+def test_resample_image_shape():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+        resample_image.main('image resampled -s 20 20 20'.split())
+        res = Image('resampled')
+
+        expv2w = transform.concat(
+            img.voxToWorldMat,
+            transform.scaleOffsetXform([0.5, 0.5, 0.5], 0))
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+        assert np.all(np.isclose(
+            np.array(transform.axisBounds(res.shape, res.voxToWorldMat)) - 0.25,
+                     transform.axisBounds(img.shape, img.voxToWorldMat)))
+
+        resample_image.main('image resampled -s 20 20 20 -o corner'.split())
+        res = Image('resampled')
+        assert np.all(np.isclose(
+            transform.axisBounds(res.shape, res.voxToWorldMat),
+            transform.axisBounds(img.shape, img.voxToWorldMat)))
+
+
+def test_resample_image_dim():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+
+        resample_image.main('image resampled -d 0.5 0.5 0.5'.split())
+
+        res = Image('resampled')
+        expv2w = transform.concat(
+            img.voxToWorldMat,
+            transform.scaleOffsetXform([0.5, 0.5, 0.5], 0))
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+
+
+def test_resample_image_ref():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+        ref = Image(make_random_image('ref.nii.gz',   dims=(20, 20, 20),
+                                      pixdims=(0.5, 0.5, 0.5)))
+
+        resample_image.main('image resampled -r ref'.split())
+
+        res    = Image('resampled')
+        expv2w = ref.voxToWorldMat
+
+        assert np.all(np.isclose(res.shape, (20, 20, 20)))
+        assert np.all(np.isclose(res.pixdim, (0.5, 0.5, 0.5)))
+        assert np.all(np.isclose(res.voxToWorldMat, expv2w))
+
+
+def test_resample_image_bad_options():
+    with tempdir():
+        img = Image(make_random_image('image.nii.gz', dims=(10, 10, 10)))
+
+        with pytest.raises(SystemExit) as e:
+            resample_image.main('image resampled -d 0.5 0.5 0.5 '
+                                '-s 20 20 20'.split())
+        assert e.value.code != 0
diff --git a/tests/test_transform.py b/tests/test_transform.py
index a6688b85f7d458765cdaf8c9ede5114403fe7cc6..70b86019aae9a7374a906dda830fd0a8f7c0966d 100644
--- a/tests/test_transform.py
+++ b/tests/test_transform.py
@@ -268,6 +268,17 @@ def test_compose_and_decompose():
     assert np.all(np.isclose(rotat, rots))
     assert np.all(np.isclose(rotaf, rmat))
 
+    # decompose should accept a 3x3
+    # affine, and return translations of 0
+    transform.decompose(xform[:3, :3])
+    sc,   of,   rot   = transform.decompose(xform[:3, :3])
+    sc,   of,   rot   = np.array(sc), np.array(of), np.array(rot)
+    assert np.all(np.isclose(sc,    [1, 1, 1]))
+    assert np.all(np.isclose(of,    [0, 0, 0]))
+    assert np.all(np.isclose(rot,   rots))
+
+
+
 
 def test_rotMatToAxisAngles(seed):