diff --git a/fsl/fslview/controls/locationpanel.py b/fsl/fslview/controls/locationpanel.py
index 91b9ab7d0dd41f988b97189faac485deabbcae6c..34fec077748d976c112af0376199b62233f605ad 100644
--- a/fsl/fslview/controls/locationpanel.py
+++ b/fsl/fslview/controls/locationpanel.py
@@ -428,30 +428,64 @@ class LocationPanel(fslpanel.FSLViewPanel):
 
         self.Freeze()
 
+
+    def _getOffsets(self, source, target, displaySpace):
+
+        # If the world location was changed, the
+        # affine->vox transformation returns voxel
+        # coordinates in the space [x - 0.5, x + 0.5].
+        # We add 0.5 to these coordinates so that the
+        # _propagate method can just floor the result
+        # to get the integer voxel coordinates. 
+        if   (source, target) == ('world', 'voxel'):
+            return (0, 0, 0), (0.5, 0.5, 0.5)
+
+        # If the reference image is being displayed in
+        # affine space, we have the same situation as
+        # above.
+        elif (source, target) == ('display', 'voxel'):
+            if displaySpace == 'affine':
+                return (0, 0, 0), (0.5, 0.5, 0.5)
+
+        # If the voxel location was changed, we want
+        # the display to be moved to the centre of the
+        # voxel
+        elif (source, target) == ('voxel', 'display'):
+            
+            if displaySpace in ('id', 'pixdim'):
+                return (0.5, 0.5, 0.5), (0, 0, 0)
+
+        return (0, 0, 0), (0, 0, 0)
+
         
     def _propagate(self, source, target, xform):
 
+        displaySpace = None
+        if self._refImage is not None:
+            opts = self._displayCtx.getOpts(self._refImage)
+            displaySpace = opts.transform
+
+        pre, post = self._getOffsets(source, target, displaySpace)
+            
         if   source == 'display': coords = self._displayCtx.location.xyz
         elif source == 'voxel':   coords = self.voxelLocation.xyz
         elif source == 'world':   coords = self.worldLocation.xyz
 
+        coords = [coords[0] + pre[0],
+                  coords[1] + pre[1],
+                  coords[2] + pre[2]]
+                
         if xform is not None: xformed = transform.transform([coords], xform)[0]
         else:                 xformed = np.array(coords)
 
+        xformed += post
+
         log.debug('Updating location ({} {} -> {} {})'.format(
             source, coords, target, xformed))
 
-        if target == 'display':
-            self._displayCtx.location.xyz = xformed
-
-        # Voxel coordinates are transformed to [x - 0.5, x + 0.5],
-        # so we floor(x + 0.5) to get the corresponding integer
-        # coordinates
-        elif target == 'voxel':
-            self.voxelLocation.xyz = np.floor(xformed + 0.5)
-            
-        elif target == 'world':
-            self.worldLocation.xyz = xformed
+        if   target == 'display': self._displayCtx.location.xyz = xformed
+        elif target == 'voxel':   self.voxelLocation.xyz = np.floor(xformed)
+        elif target == 'world':   self.worldLocation.xyz = xformed
         
     
     def _postPropagate(self):
@@ -532,11 +566,15 @@ class LocationPanel(fslpanel.FSLViewPanel):
                     [self._displayCtx.location.xyz],
                     opts.getTransform('display', 'voxel'))[0]
 
-                # The above transformation gives us
+                # When displaying in world/affine space,
+                # the above transformation gives us
                 # values between [x - 0.5, x + 0.5] for
                 # voxel x, so we need to floor(x + 0.5)
                 # to get the actual voxel coordinates
-                vloc = tuple(map(int, np.floor(vloc + 0.5)))
+                if opts.transform == 'affine': vloc = np.floor(vloc + 0.5)
+                else:                          vloc = np.floor(vloc)
+
+                vloc = tuple(map(int, vloc))
 
                 if overlay.is4DImage():
                     vloc = vloc + (opts.volume,)
diff --git a/fsl/fslview/displaycontext/volumeopts.py b/fsl/fslview/displaycontext/volumeopts.py
index 0046cb69bd9221c097f55c0ac44fe4ebe9a75df6..9495186b093d3841acd79e427dfde826690ea105 100644
--- a/fsl/fslview/displaycontext/volumeopts.py
+++ b/fsl/fslview/displaycontext/volumeopts.py
@@ -34,7 +34,8 @@ class ImageOpts(fsldisplay.DisplayOpts):
 
     
     resolution = props.Real(maxval=10, default=1, clamped=True)
-    """Data resolution in world space. The minimum value is set in __init__.""" 
+    """Data resolution in world space. The minimum value is set in __init__.
+    """ 
 
 
     transform = props.Choice(
@@ -43,8 +44,12 @@ class ImageOpts(fsldisplay.DisplayOpts):
                 strings.choices['ImageOpts.transform.pixdim'],
                 strings.choices['ImageOpts.transform.id']],
         default='pixdim')
-    """This property defines how the overlay should be transformd into the display
-    coordinate system.
+    """This property defines how the overlay should be transformd into the
+    display coordinate system.
+
+    note:: Write a note here about voxels in ``affine`` mapping to space
+           ``[x - 0.5, x + 0.5]``, whereas the other two mapping to
+           ``[x, x + 1]``.
     
       - ``affine``: Use the affine transformation matrix stored in the image
         (the ``qform``/``sform`` fields in NIFTI1 headers).
@@ -100,8 +105,13 @@ class ImageOpts(fsldisplay.DisplayOpts):
         :attr:`.DisplayOpts.bounds` property accordingly.
         """
 
+        if self.transform == 'affine': origin = 'centre'
+        else:                          origin = 'corner'
+
         lo, hi = transform.axisBounds(
-            self.overlay.shape[:3], self.getTransform('voxel', 'display'))
+            self.overlay.shape[:3],
+            self.getTransform('voxel', 'display'),
+            origin=origin)
         
         self.bounds[:] = [lo[0], hi[0], lo[1], hi[1], lo[2], hi[2]]
 
diff --git a/fsl/fslview/gl/annotations.py b/fsl/fslview/gl/annotations.py
index 383ba7852eb2197872a72d74f1809ea0b2a99051..5e39f871c1ef92dee116cc3c782f8e8852a6d731 100644
--- a/fsl/fslview/gl/annotations.py
+++ b/fsl/fslview/gl/annotations.py
@@ -438,7 +438,8 @@ class VoxelSelection(AnnotationObject):
                                             yax,
                                             zpos,
                                             self.voxToDisplayMat,
-                                            self.displayToVoxMat)
+                                            self.displayToVoxMat,
+                                            origin='corner')
 
         verts = np.array(verts, dtype=np.float32).ravel('C')
         texs  = np.array(texs,  dtype=np.float32).ravel('C')
diff --git a/fsl/fslview/gl/globject.py b/fsl/fslview/gl/globject.py
index 611118e1ff76c3353c80c647826cc7ae785f5a30..39a5fe904d6c17b61d8de8b278efbe1d6975f87a 100644
--- a/fsl/fslview/gl/globject.py
+++ b/fsl/fslview/gl/globject.py
@@ -294,14 +294,19 @@ class GLImageObject(GLObject):
         
           - A ``6*3 numpy.float32`` array containing the texture coordinates
             corresponding to each vertex
-        """ 
+        """
+
+        if self.displayOpts.transform == 'affine': origin = 'centre'
+        else:                                      origin = 'corner'
+
         vertices, voxCoords, texCoords = glroutines.slice2D(
             self.image.shape[:3],
             self.xax,
             self.yax,
             zpos, 
             self.displayOpts.getTransform('voxel',   'display'),
-            self.displayOpts.getTransform('display', 'voxel'))
+            self.displayOpts.getTransform('display', 'voxel'),
+            origin=origin)
 
         if xform is not None: 
             vertices = transform.transform(vertices, xform)
diff --git a/fsl/fslview/gl/routines.py b/fsl/fslview/gl/routines.py
index f47810ead95f3233944984885e6c62ae77db2fb8..3cba3576b14dbefff5a4ffae98cc8bb380f9e8c2 100644
--- a/fsl/fslview/gl/routines.py
+++ b/fsl/fslview/gl/routines.py
@@ -347,7 +347,8 @@ def slice2D(dataShape,
             zpos,
             voxToDisplayMat,
             displayToVoxMat,
-            geometry='triangles'):
+            geometry='triangles',
+            origin='centre'):
     """Generates and returns vertices which denote a slice through an
     array of the given ``dataShape``, parallel to the plane defined by the
     given ``xax`` and ``yax`` and at the given z position, in the space
@@ -373,6 +374,13 @@ def slice2D(dataShape,
         |   |
         0---1
 
+    If ``origin`` is set to ``centre`` (the default), it is assumed that
+    a voxel at location ``(x, y, z)`` is located in the space
+    ``(x - 0.5 : x + 0.5, y - 0.5 : y + 0.5, z - 0.5 : z + 0.5). Otherwise,
+    if ``origin`` is set to ``corner``, a voxel at location ``(x, y, z)``
+    is assumed to be located in the space
+    ``(x : x + 1, y : y + 1, z : z + 1)``.
+
     :arg dataShape:       Number of elements along each dimension in the
                           image data.
     
@@ -389,6 +397,11 @@ def slice2D(dataShape,
                           system.
 
     :arg displayToVoxMat: Inverse of the ``voxToDisplayMat``.
+
+    :arg geometry:        ``square`` or ``triangle``.
+
+    :arg origin:          ``centre`` or ``corner``. See the
+                          :func:`.transform.axisBounds` function.
     
     Returns a tuple containing:
     
@@ -401,12 +414,11 @@ def slice2D(dataShape,
 
       - A ``N*3`` ``numpy.float32`` array containing the texture coordinates
         that correspond to the voxel coordinates.
-
     """
 
     zax        = 3 - xax - yax
-    xmin, xmax = transform.axisBounds(dataShape, voxToDisplayMat, xax)
-    ymin, ymax = transform.axisBounds(dataShape, voxToDisplayMat, yax)
+    xmin, xmax = transform.axisBounds(dataShape, voxToDisplayMat, xax, origin)
+    ymin, ymax = transform.axisBounds(dataShape, voxToDisplayMat, yax, origin)
 
     if geometry == 'triangles':
 
@@ -437,7 +449,9 @@ def slice2D(dataShape,
     # default centered at 0 (i.e. the space of a voxel
     # lies in the range [-0.5, 0.5]), but we want voxel
     # coordinates to map to the effective range [0, 1]
-    voxCoords = voxCoords + 0.5
+    if origin == 'centre':
+        voxCoords = voxCoords + 0.5
+        
     texCoords = voxCoords / dataShape
 
     return vertices, voxCoords, texCoords
diff --git a/fsl/fslview/profiles/orthoeditprofile.py b/fsl/fslview/profiles/orthoeditprofile.py
index 8e366f3d98e28898fc1d2878e5293e015938a29e..03864948ac3760b196cc8cb0d05e988c46182887 100644
--- a/fsl/fslview/profiles/orthoeditprofile.py
+++ b/fsl/fslview/profiles/orthoeditprofile.py
@@ -303,10 +303,7 @@ class OrthoEditProfile(orthoviewprofile.OrthoViewProfile):
         voxel = transform.transform(
             [canvasPos], opts.getTransform('display', 'voxel'))[0]
 
-        # Using floor(voxel+0.5) because, when at the
-        # midpoint, I want to round up. np.round rounds
-        # to the nearest even number, which is not ideal
-        voxel = np.array(np.floor(voxel + 0.5), dtype=np.int32)
+        voxel = np.int32(np.floor(voxel))
 
         return voxel
 
@@ -337,12 +334,13 @@ class OrthoEditProfile(orthoviewprofile.OrthoViewProfile):
         if blockSize is None:
             blockSize = self.selectionSize
 
-        lo = [(v - 0.5) - int(np.floor((blockSize - 1) / 2.0)) for v in voxel]
-        hi = [(v + 0.5) + int(np.ceil(( blockSize - 1) / 2.0)) for v in voxel]
+        # TODO Double check this, as it seems to be wrong
+        lo = [(v)     - int(np.floor((blockSize - 1) / 2.0)) for v in voxel]
+        hi = [(v + 1) + int(np.ceil(( blockSize - 1) / 2.0)) for v in voxel]
 
         if not self.selectionIs3D:
-            lo[canvas.zax] = voxel[canvas.zax] - 0.5
-            hi[canvas.zax] = voxel[canvas.zax] + 0.5
+            lo[canvas.zax] = voxel[canvas.zax]
+            hi[canvas.zax] = voxel[canvas.zax] + 1
 
         corners       = np.zeros((8, 3))
         corners[0, :] = lo[0], lo[1], lo[2]
diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py
index b77526f25ae9a5704fdf7c59996d9d15b6323ed8..2831152408127552c46bf79715938271de493905 100644
--- a/fsl/utils/transform.py
+++ b/fsl/utils/transform.py
@@ -50,10 +50,11 @@ def scaleOffsetXform(scales, offsets):
     return xform
 
 
-def axisBounds(shape, xform, axes=None):
+def axisBounds(shape, xform, axes=None, origin='centre'):
     """Returns the (lo, hi) bounds of the specified axis/axes.
-
-    This function assumes that voxel indices correspond to the voxel
+    
+    If the ``origin`` parameter is set to  ``centre`` (the default),
+    this function assumes that voxel indices correspond to the voxel
     centre. For example, the voxel at ``(4, 5, 6)`` covers the space:
     
       ``[3.5 - 4.5, 4.5 - 5.5, 5.5 - 6.5]``
@@ -65,8 +66,25 @@ def axisBounds(shape, xform, axes=None):
     to the corner at
 
     ``(shape[0] - 0.5, shape[1] - 0.5, shape[1] - 0.5)``
+
+    If the ``origin`` parameter is set to ``corner``, this function
+    assumes that voxel indices correspond to the voxel corner. In this
+    case, a voxel  at ``(4, 5, 6)`` covers the space:
+    
+      ``[4 - 5, 5 - 6, 6 - 7]``
+    
+    So the bounds of the specified shape extends from the corner at
+
+    ``(0, 0, 0)``
+
+    to the corner at
+
+    ``(shape[0], shape[1]5, shape[1])``.
     """
 
+    if origin not in ('centre', 'corner'):
+        raise ValueError('Invalid origin value: {}'.format(origin))
+
     scalar = False
 
     if axes is None:
@@ -80,18 +98,26 @@ def axisBounds(shape, xform, axes=None):
 
     points = np.zeros((8, 3), dtype=np.float32)
 
-    x -= 0.5
-    y -= 0.5
-    z -= 0.5
-
-    points[0, :] = [-0.5, -0.5, -0.5]
-    points[1, :] = [-0.5, -0.5,  z]
-    points[2, :] = [-0.5,  y,   -0.5]
-    points[3, :] = [-0.5,  y,    z]
-    points[4, :] = [x,    -0.5, -0.5]
-    points[5, :] = [x,    -0.5,  z]
-    points[6, :] = [x,     y,   -0.5]
-    points[7, :] = [x,     y,    z]
+    if origin == 'centre':
+        x0 = -0.5
+        y0 = -0.5
+        z0 = -0.5
+        x -=  0.5
+        y -=  0.5
+        z -=  0.5
+    else:
+        x0 = 0
+        y0 = 0
+        z0 = 0
+
+    points[0, :] = [x0, y0, z0]
+    points[1, :] = [x0, y0,  z]
+    points[2, :] = [x0,  y, z0]
+    points[3, :] = [x0,  y,  z]
+    points[4, :] = [x,  y0, z0]
+    points[5, :] = [x,  y0,  z]
+    points[6, :] = [x,   y, z0]
+    points[7, :] = [x,   y,  z]
 
     tx = transform(points, xform)