diff --git a/fsl/data/image.py b/fsl/data/image.py
index 6c55b2e6c470fb83bd641c4d94d94900a8f91b5c..fe179c4e714f5a2c6bfbd2cf609bcd573862ae7c 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -225,6 +225,20 @@ class Nifti1(object):
         return int(code)
 
 
+    def axisMapping(self, xform):
+        """Returns the (approximate) correspondence of each axis in the source
+        coordinate system to the axes in the destination coordinate system,
+        where the source and destinations are defined by the given affine
+        transformation matrix.
+        """
+
+        import nibabel as nib
+
+        inaxes = [[-1, 1], [-2, 2], [-3, 3]]
+
+        return nib.orientations.aff2axcodes(xform, inaxes)
+
+
     def getOrientation(self, axis, xform):
         """Returns a code representing the orientation of the specified data
         axis in the coordinate system defined by the given transformation
diff --git a/fsl/fsleyes/displaycontext/volumeopts.py b/fsl/fsleyes/displaycontext/volumeopts.py
index 9a4ac675118d451773e1f1f48866a3c4c5abdf49..e26646a3a95a8eef0b10cbf3df17ba71ceb3d82f 100644
--- a/fsl/fsleyes/displaycontext/volumeopts.py
+++ b/fsl/fsleyes/displaycontext/volumeopts.py
@@ -369,7 +369,7 @@ class Nifti1Opts(fsldisplay.DisplayOpts):
         return (0, 0, 0), (0, 0, 0)
 
 
-    def transformCoords(self, coords, from_, to_):
+    def transformCoords(self, coords, from_, to_, vround=False):
         """Transforms the given coordinates from ``from_`` to ``to_``.
 
         The ``from_`` and ``to_`` parameters must both be one of:
@@ -379,41 +379,85 @@ class Nifti1Opts(fsldisplay.DisplayOpts):
            - ``world``:   The image world coordinate system
            - ``custom``:  The coordinate system defined by the custom
                           transformation matrix (see :attr:`customXform`)
+
+        :arg coords: Coordinates to transform
+        :arg from_:  Space to transform from
+        :arg to_:    Space to transform to
+        :arg vround: If ``True``, and ``to_ == 'voxel'``, the transformed
+                     coordinates are rounded to the nearest integer.
         """
 
         xform     = self.getTransform(       from_, to_)
         pre, post = self.getTransformOffsets(from_, to_)
         
-        coords    = np.array(coords) + pre
-        coords    = transform.transform(coords, xform)
-
-        return coords + post
+        coords    = np.array(coords)                   + pre
+        coords    = transform.transform(coords, xform) + post
+        
+        # Round to integer voxel coordinates?
+        if to_ == 'voxel' and vround:
+
+            # The transformation matrices treat integer
+            # voxel coordinates as the voxel centre (e.g.
+            # a voxel [3, 4, 5] fills the space:
+            # 
+            # [2.5-3.5, 3.5-4.5, 4.5-5.5].
+            #
+            # So all we need to do is round to the
+            # nearest integer.
+            #
+            # Note that the numpy.round function breaks
+            # ties (e.g. 7.5) by rounding to the nearest
+            # *even* integer, which can cause funky
+            # behaviour. So we take (floor(x)+0.5) instead
+            # of rounding, to force consistent behaviour
+            # (i.e. always rounding central values up).
+            coords = np.floor(coords + 0.5)
+
+        return coords
 
     
-    def getVoxel(self):
+    def getVoxel(self, xyz=None, clip=True, vround=True):
         """Calculates and returns the voxel coordinates corresponding to the
-        current :attr:`.DisplayContext.location` for the :class:`.Nifti1`
-        associated with this ``Nifti1Opts`` instance.
+        given location (assumed to be in the display coordinate system) for
+        the :class:`.Nifti1` associated with this ``Nifti1Opts`` instance..
+
+        :arg xyz:    Display space location to convert to voxels. If not
+                     provided, the current :attr:`.DisplayContext.location`
+                     is used.
+
+        :arg clip:   If ``False``, and the transformed coordinates are out of
+                     the voxel coordinate bounds, the coordinates returned
+                     anyway. Defaults to ``True``.
 
-        Returns ``None`` if the current location is outside of the image
-        bounds.
+
+        :arg vround: If ``True``, the returned voxel coordinates are rounded
+                     to the nearest integer. Otherwise they may be fractional.
+                    
+
+        :returns:    ``None`` if the location is outside of the image bounds,
+                     unless ``clip=False``.
         """
 
+        if xyz is not None: x, y, z = xyz
+        else:               x, y, z = self.displayCtx.location.xyz
+
         overlay = self.overlay
-        x, y, z = self.displayCtx.location.xyz
+        vox     = self.transformCoords([[x, y, z]],
+                                       'display',
+                                       'voxel',
+                                       vround=vround)[0]
+
+        if vround:
+            vox = map(int, vox)
 
-        vox     = self.transformCoords([[x, y, z]], 'display', 'voxel')[0]
-        vox     = map(int, np.round(vox))
+        if not clip:
+            return vox
 
-        if vox[0] < 0                 or \
-           vox[1] < 0                 or \
-           vox[2] < 0                 or \
-           vox[0] >= overlay.shape[0] or \
-           vox[1] >= overlay.shape[1] or \
-           vox[2] >= overlay.shape[2]:
-            return None
+        for ax in (0, 1, 2):
+            if vox[ax] < 0 or vox[ax] >= overlay.shape[ax]:
+                return None
 
-        return vox 
+        return vox
 
 
     def displayToStandardCoordinates(self, coords):
diff --git a/fsl/fsleyes/gl/routines.py b/fsl/fsleyes/gl/routines.py
index f3d7ec1b3755ec8e48278f29ed912ec2d6c6cbdd..110917848fee656c5cc5b319c002bcd800a16286 100644
--- a/fsl/fsleyes/gl/routines.py
+++ b/fsl/fsleyes/gl/routines.py
@@ -489,8 +489,17 @@ def slice2D(dataShape,
     # coordinates to map to the effective range [0, 1]
     if origin == 'centre':
         voxCoords = voxCoords + 0.5
+        texCoords = voxCoords / dataShape
+
+    # If the origin is the voxel corner, we still need
+    # to offset the texture coordinates, as the voxel
+    # coordinate range [-0.5, shape-0.5] must be scaled
+    # to the texture coordinate range [0.0, 1.0].
+    elif origin == 'corner':
+        texCoords = (voxCoords + 0.5) / dataShape
         
-    texCoords = voxCoords / dataShape
+    else:
+        raise ValueError('Unrecognised origin: {}'.format(origin))
 
     return vertices, voxCoords, texCoords
 
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index 095b645cf4fa0a23da317949c04281f986ee4b53..ae5b8b22f18bbd677e53c01a208862cf9503acab 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -70,9 +70,14 @@ def scaleOffsetXform(scales, offsets):
     return xform
 
 
-def axisBounds(shape, xform, axes=None, origin='centre'):
-    """Returns the ``(lo, hi)`` bounds of the specified axis/axes in the world
-    coordinate system defined by ``xform``.
+def axisBounds(shape,
+               xform,
+               axes=None,
+               origin='centre',
+               boundary='high',
+               offset=1e-4):
+    """Returns the ``(lo, hi)`` bounds of the specified axis/axes in the
+    world coordinate system defined by ``xform``.
     
     If the ``origin`` parameter is set to  ``centre`` (the default),
     this function assumes that voxel indices correspond to the voxel
@@ -100,21 +105,35 @@ def axisBounds(shape, xform, axes=None, origin='centre'):
 
     to the corner at
 
-      ``(shape[0], shape[1]5, shape[1])``.
+      ``(shape[0], shape[1], shape[1])``.
 
+
+    If the ``boundary`` parameter is set to ``high``, the high voxel bounds
+    are reduced by a small amount (specified by the ``offset`` parameter)
+    before they are transformed to the world coordinate system.  If
+    ``boundary`` is set to ``low``, the low bounds are increased by a small
+    amount.  The ``boundary`` parameter can also be set to ``'both'``, or
+    ``None``. This option is provided so that you can ensure that the
+    resulting bounds will always be contained within the image space.
     
-    :arg shape:  The ``(x, y, z)`` shape of the data.
+    :arg shape:    The ``(x, y, z)`` shape of the data.
+
+    :arg xform:    Transformation matrix which transforms voxel coordinates
+                   to the world coordinate system.
+
+    :arg axes:     The world coordinate system axis bounds to calculate.
 
-    :arg xform:  Transformation matrix which transforms voxel coordinates
-                 to the world coordinate system.
+    :arg origin:   Either ``'centre'`` (the default) or ``'origin'``.
 
-    :arg axes:   The world coordinate system axis bounds to calculate.
+    :arg boundary: Either ``'high'`` (the default), ``'low'``, ''`both'``,
+                   or ``None``. 
 
-    :arg origin: Either ``'centre'`` or ``'origin'``
+    :arg offset:   Amount by which the boundary voxel coordinates should be
+                   offset. Defaults to ``1e-4``.
 
-    :returns:    A list of tuples, one for each axis specified in the ``axes``
-                 argument. Each tuple contains the ``(lo, hi)`` bounds of the
-                 corresponding world coordinate system axis.
+    :returns:      A list of tuples, one for each axis specified in the 
+                   ``axes`` argument. Each tuple contains the ``(lo, hi)`` 
+                   bounds of the corresponding world coordinate system axis.
     """
 
     origin = origin.lower()
@@ -151,6 +170,16 @@ def axisBounds(shape, xform, axes=None, origin='centre'):
         y0 = 0
         z0 = 0
 
+    if boundary in ('low', 'both'):
+        x0 += offset
+        y0 += offset
+        z0 += offset
+        
+    if boundary in ('high', 'both'):
+        x  -= offset
+        y  -= offset
+        z  -= offset
+
     points[0, :] = [x0, y0, z0]
     points[1, :] = [x0, y0,  z]
     points[2, :] = [x0,  y, z0]
@@ -168,7 +197,6 @@ def axisBounds(shape, xform, axes=None, origin='centre'):
     if scalar: return (lo[0], hi[0])
     else:      return (lo,    hi)
 
-
         
 def transform(p, xform, axes=None):
     """Transforms the given set of points ``p`` according to the given affine