From 0da23afc26632cc832a4e141d3cd6c6af4817ad8 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauld.mccarthy@gmail.com>
Date: Wed, 10 Feb 2016 11:25:04 +0000
Subject: [PATCH] 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                        | 14 ++++
 fsl/fsleyes/displaycontext/volumeopts.py | 86 ++++++++++++++++++------
 fsl/fsleyes/gl/routines.py               | 11 ++-
 fsl/utils/transform.py                   | 54 +++++++++++----
 4 files changed, 130 insertions(+), 35 deletions(-)

diff --git a/fsl/data/image.py b/fsl/data/image.py
index 6c55b2e6c..fe179c4e7 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 9a4ac6751..e26646a3a 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 f3d7ec1b3..110917848 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 095b645cf..ae5b8b22f 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
-- 
GitLab