VolumeOpts.transformCoords will optionally round coordinates (if the

target space is 'voxel'), and will optionally not clip the
result. VolumeOpts.getVoxel allows coordinates to be passed
in. Transform.axisBounds allows an offset to be applied to voxel
coordinates. New method, Image.axisMapping, which basically wraps
nib.aff2axcodes.

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

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):

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
 

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