From d9fca5c309908795b52f55ebef2282e07c87c334 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauld.mccarthy@gmail.com>
Date: Thu, 13 Aug 2015 17:53:24 +0100
Subject: [PATCH] Gaargh, I've confused myself badly. Images in id/pixdim space
should be displayed with the bottom left voxel at display location (0, 0).
This means that calculating voxel location from display coordinates is
different depending on the space that the image is displayed in. I think I've
figured out the location panel, but need to double check, and need to take a
closer look at the editor.
---
fsl/fslview/controls/locationpanel.py | 64 +++++++++++++++++++-----
fsl/fslview/displaycontext/volumeopts.py | 18 +++++--
fsl/fslview/gl/annotations.py | 3 +-
fsl/fslview/gl/globject.py | 9 +++-
fsl/fslview/gl/routines.py | 24 +++++++--
fsl/fslview/profiles/orthoeditprofile.py | 14 +++---
fsl/utils/transform.py | 56 +++++++++++++++------
7 files changed, 140 insertions(+), 48 deletions(-)
diff --git a/fsl/fslview/controls/locationpanel.py b/fsl/fslview/controls/locationpanel.py
index 91b9ab7d0..34fec0777 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 0046cb69b..9495186b0 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 383ba7852..5e39f871c 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 611118e1f..39a5fe904 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 f47810ead..3cba3576b 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 8e366f3d9..03864948a 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 b77526f25..283115240 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)
--
GitLab