From d65afb36fe2ace31d5e16b2b78af12cfdbbfa200 Mon Sep 17 00:00:00 2001 From: Paul McCarthy <pauld.mccarthy@gmail.com> Date: Wed, 22 Oct 2014 14:14:37 +0100 Subject: [PATCH] All the code is broken. Am in process of decoupling display and world coordinate systemss, with the ultimate aim of having the location panel always reporting image world coordinates, regardless of the way in which images are being displayed. --- fsl/data/image.py | 286 ++++------------------------------ fsl/fslview/displaycontext.py | 171 ++++++++++++++++---- fsl/utils/transform.py | 118 ++++++++++++++ 3 files changed, 284 insertions(+), 291 deletions(-) create mode 100644 fsl/utils/transform.py diff --git a/fsl/data/image.py b/fsl/data/image.py index 00bfd8e4f..42fa6944f 100644 --- a/fsl/data/image.py +++ b/fsl/data/image.py @@ -12,16 +12,14 @@ import sys import logging import tempfile import collections -import subprocess as sp -import os.path as op +import subprocess as sp +import os.path as op -import numpy as np -import numpy.linalg as linalg -import nibabel as nib +import numpy as np +import nibabel as nib import props - -import fsl.data.imagefile as imagefile +import fsl.data.imagefile as imagefile log = logging.getLogger(__name__) @@ -48,6 +46,7 @@ NIFTI_XFORM_MNI_152 = 4 # image is being displayed in voxel space NIFTI_XFORM_VOXEL = 5 + def _loadImageFile(filename): """Given the name of an image file, loads it using nibabel. @@ -120,30 +119,26 @@ class Image(props.HasProperties): The following attributes are present on an :class:`Image` object: - :ivar nibImage: The :mod:`nibabel` image object. - - :ivar data: A reference to the image data, stored as a - :mod`numpy` array. + :ivar nibImage: The :mod:`nibabel` image object. - :ivar shape: A list/tuple containing the number of voxels - along each image dimension. + :ivar data: A reference to the image data, stored as a + :mod`numpy` array. - :ivar pixdim: A list/tuple containing the size of one voxel - along each image dimension. + :ivar shape: A list/tuple containing the number of voxels + along each image dimension. - :ivar voxToWorldMat: A 4*4 array specifying the affine transformation - for transforming voxel coordinates into real world - coordinates. + :ivar pixdim: A list/tuple containing the size of one voxel + along each image dimension. - :ivar worldToVoxMat: A 4*4 array specifying the affine transformation - for transforming real world coordinates into voxel - coordinates. - - :ivar imageFile: The name of the file that the image was loaded from. + :ivar voxToWorldMat: A 4*4 array specifying the affine transformation + for transforming voxel coordinates into real world + coordinates. + + :ivar imageFile: The name of the file that the image was loaded from. - :ivar tempFile: The name of the temporary file which was created (in - the event that the image was large and was gzipped - - see :func:`_loadImageFile`). + :ivar tempFile: The name of the temporary file which was created (in + the event that the image was large and was gzipped - + see :func:`_loadImageFile`). """ @@ -156,27 +151,6 @@ class Image(props.HasProperties): """This property defines the type of image data.""" - transform = props.Choice( - collections.OrderedDict([ - ('affine', 'Use qform/sform transformation matrix'), - ('pixdim', 'Use pixdims only'), - ('id', 'Do not use qform/sform or pixdims')]), - default='pixdim') - """This property defines how the image should be transformd into real world - space. - - - ``affine``: Use the affine transformation matrix stored in the image - (the ``qform``/``sform`` fields in NIFTI1 headers). - - - ``pixdim``: Scale voxel sizes by the ``pixdim`` fields in the image - header. - - - ``id``: Perform no scaling or transformation - voxels will be - interpreted as :math:`1mm^3` isotropic, with the origin at voxel - (0,0,0). - """ - - name = props.String() """The name of this image.""" @@ -221,16 +195,10 @@ class Image(props.HasProperties): self.tempFile = None self.imageFile = None - self.data = self.nibImage.get_data() - self.shape = self.nibImage.get_shape() - self.pixdim = self.nibImage.get_header().get_zooms() - - self.addListener( - 'transform', - '{}_{}'.format(self.__class__.__name__, self.name), - self._transformChanged) - - self._transformChanged() + self.data = self.nibImage.get_data() + self.shape = self.nibImage.get_shape() + self.pixdim = self.nibImage.get_header().get_zooms() + self.voxToWorldMat = np.array(self.nibImage.get_affine()) if len(self.shape) < 3 or len(self.shape) > 4: raise RuntimeError('Only 3D or 4D images are supported') @@ -253,140 +221,7 @@ class Image(props.HasProperties): def is4DImage(self): """Returns ``True`` if this image is 4D, ``False`` otherwise. """ - return len(self.shape) > 3 and self.shape[3] > 1 - - - def _transformChanged(self, *a): - """This method is called when the :attr:`transform` property value - changes. It updates the :attr:`voxToWorldMat`, :attr:`worldToVoxMat`, - and :attr:`pixdim` attributes to reflect the new transformation - type. - """ - - if self.transform == 'affine': - voxToWorldMat = self.nibImage.get_affine() - elif self.transform == 'pixdim': - pixdim = self.nibImage.get_header().get_zooms() - voxToWorldMat = np.diag([pixdim[0], pixdim[1], pixdim[2], 1.0]) - elif self.transform == 'id': - voxToWorldMat = np.identity(4) - - self.voxToWorldMat = np.array(voxToWorldMat, dtype=np.float32) - self.worldToVoxMat = linalg.inv(self.voxToWorldMat) - - self.voxToWorldMat = self.voxToWorldMat.transpose() - self.worldToVoxMat = self.worldToVoxMat.transpose() - - if self.transform == 'affine': - pixdim = [self.axisLength(ax) / self.shape[ax] for ax in range(3)] - elif self.transform == 'pixdim': - pixdim = self.nibImage.get_header().get_zooms() - elif self.transform == 'id': - pixdim = [1.0, 1.0, 1.0] - - self.pixdim = pixdim - - # for pixdim/identity transformations, we want the world - # location (0, 0, 0) to map to voxel location (0, 0, 0) - if self.transform in ['pixdim', 'id']: - for i in range(3): - self.voxToWorldMat[3, i] = self.pixdim[i] * 0.5 - self.worldToVoxMat[3, i] = -0.5 - - log.debug('Image {} transformation matrix changed: {}'.format( - self.name, self.voxToWorldMat)) - log.debug('Inverted matrix: {}'.format(self.worldToVoxMat)) - - - def imageBounds(self, axis): - """Return the bounds (min, max) of the image, in real world - - - The returned bounds give the coordinates, along the specified axis, of - a bounding box which contains the entire image. - """ - - x, y, z = self.shape[:3] - - x -= 0.5 - y -= 0.5 - z -= 0.5 - - points = np.zeros((8, 3), dtype=np.float32) - - 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] - - - tx = self.voxToWorld(points) - - lo = tx[:, axis].min() - hi = tx[:, axis].max() - - return (lo, hi) - - - def axisLength(self, axis): - """Return the length, in real world units, of the specified axis. - """ - - points = np.zeros((2, 3), dtype=np.float32) - points[:] = [-0.5, -0.5, -0.5] - points[1, axis] = self.shape[axis] - 0.5 - - tx = self.voxToWorld(points) - - # euclidean distance between each boundary point - return sum((tx[0, :] - tx[1, :]) ** 2) ** 0.5 - - - def worldToVox(self, p, axes=None): - """Transforms the given set of points in voxel coordinates to points - in world coordinates, according to the current :attr:`transform`. - - The returned array is either a :class:`numpy.float64` array, or a - single ``float`` value, depending on the input. There is no guarantee - that the returned array of voxel coordinates is within the bounds of - the image shape. Parameters: - - :arg p: N*A array, where N is the number of points, and A - is the number of axes to consider (default: 3). - - :arg axes: If ``None``, it is assumed that the input p is a N*3 - array, with each point being specified by x,y,z - coordinates. If a single value in the range (0-2), - it is assumed that p is a 1D array. Or, if a - sequence of 2 or 3 values, p must be an array of - N*2 or N*3, respectively. - """ - - voxp = self._transform(p, self.worldToVoxMat, axes) - voxp = np.array(voxp, dtype=np.float64) - - if voxp.size == 1: return voxp[0] - else: return voxp - - - def voxToWorld(self, p, axes=None): - """Transforms the given set of points in world coordinates to - points in voxel coordinates, according to the current - :attr:`transform`. - - The returned array is either a :class:`numpy.float64` array, - or a single ``float`` value, depending on the input. See the - :meth:`worldToVox` method for more details. - """ - - worldp = self._transform(p, self.voxToWorldMat, axes) - - if worldp.size == 1: return float(worldp) - else: return worldp + return len(self.shape) > 3 and self.shape[3] > 1 def getXFormCode(self): @@ -459,64 +294,7 @@ class Image(props.HasProperties): (ORIENT_S2I, ORIENT_I2S)))[axis] return code - - def _transform(self, p, a, axes): - """Used by the :meth:`worldToVox` and :meth:`voxToWorld` methods. - - Transforms the given set of points ``p`` according to the given affine - transformation ``a``. The transformed points are returned as a - :class:``numpy.float64`` array. - """ - - p = self._fillPoints(p, axes) - t = np.zeros((len(p), 3), dtype=np.float64) - - x = p[:, 0] - y = p[:, 1] - z = p[:, 2] - - t[:, 0] = x * a[0, 0] + y * a[1, 0] + z * a[2, 0] + a[3, 0] - t[:, 1] = x * a[0, 1] + y * a[1, 1] + z * a[2, 1] + a[3, 1] - t[:, 2] = x * a[0, 2] + y * a[1, 2] + z * a[2, 2] + a[3, 2] - - if axes is None: axes = [0, 1, 2] - - return t[:, axes] - - - def _fillPoints(self, p, axes): - """Used by the :meth:`_transform` method. Turns the given array p into - a N*3 array of x,y,z coordinates. The array p may be a 1D array, or an - N*2 or N*3 array. - """ - - if not isinstance(p, collections.Iterable): p = [p] - - p = np.array(p) - - if axes is None: return p - - if not isinstance(axes, collections.Iterable): axes = [axes] - - if p.ndim == 1: - p = p.reshape((len(p), 1)) - - if p.ndim != 2: - raise ValueError('Points array must be either one or two ' - 'dimensions') - - if len(axes) != p.shape[1]: - raise ValueError('Points array shape does not match specified ' - 'number of axes') - - newp = np.zeros((len(p), 3), dtype=p.dtype) - - for i, ax in enumerate(axes): - newp[:, ax] = p[:, i] - - return newp - - + def getAttribute(self, name): """Retrieve the attribute with the given name. @@ -564,16 +342,6 @@ class ImageList(props.HasProperties): """A list of :class:`Image` objects. to be displayed""" - bounds = props.Bounds(ndims=3) - """This property contains the min/max values of - a bounding box (in real world coordinates) which - is big enough to contain all of the images in the - :attr:`images` list. This property shouid be - read-only, but I don't have a way to enforce it - (yet). - """ - - def __init__(self, images=None): """Create an ImageList object from the given sequence of :class:`Image` objects.""" diff --git a/fsl/fslview/displaycontext.py b/fsl/fslview/displaycontext.py index 60af89a5e..118a9ddb4 100644 --- a/fsl/fslview/displaycontext.py +++ b/fsl/fslview/displaycontext.py @@ -5,13 +5,16 @@ # Author: Paul McCarthy <pauldmccarthy@gmail.com> # -from collections import OrderedDict +import sys +import collections +from collections import OrderedDict import numpy as np import props -import fsl.data.image as fslimage +import fsl.data.image as fslimage +import fsl.utils.transform as transform import fsl.fslview.colourmaps as fslcm class ImageDisplay(props.HasProperties): @@ -76,11 +79,25 @@ class ImageDisplay(props.HasProperties): volume = props.Int(minval=0, maxval=0, default=0, clamped=True) """If a 4D image, the current volume to display.""" - - transform = fslimage.Image.transform - """How the image is transformed from voxel to real world coordinates. - This property is bound to the :attr:`~fsl.data.image.Image.transform` - property of the image associated with this :class:`ImageDisplay`. + + transform = props.Choice( + collections.OrderedDict([ + ('affine', 'Use qform/sform transformation matrix'), + ('pixdim', 'Use pixdims only'), + ('id', 'Do not use qform/sform or pixdims')]), + default='pixdim') + """This property defines how the image should be transformd into the display + coordinate system. + + - ``affine``: Use the affine transformation matrix stored in the image + (the ``qform``/``sform`` fields in NIFTI1 headers). + + - ``pixdim``: Scale voxel sizes by the ``pixdim`` fields in the image + header. + + - ``id``: Perform no scaling or transformation - voxels will be + interpreted as :math:`1mm^3` isotropic, with the origin at voxel + (0,0,0). """ @@ -241,11 +258,9 @@ class ImageDisplay(props.HasProperties): self.image = image - # bind self.transform and self.name to - # image.transform/image.name, so changes + # bind self.name to # image.name, so changes # in one are propagated to the other - self.bindProps('transform', image) - self.bindProps('name', image) + self.bindProps('name', image) # Attributes controlling image display. Only # determine the real min/max for small images - @@ -275,15 +290,49 @@ class ImageDisplay(props.HasProperties): if image.is4DImage(): self.setConstraint('volume', 'maxval', image.shape[3] - 1) + self.addListener( + 'transform', + 'ImageDisplay_{}'.format(id(self)), + self._transformChanged) + + self._transformChanged() + + + def _transformChanged(self, *a): + """Called when the :attr:`transform` property is changed. + + Generates transformation matrices for transforming between voxel and + display coordinate space. + + If :attr:`transform` is set to ``affine``, the :attr:`interpolation` + property is changed to ``spline. Otherwise, it is set to ``none``. + """ + + if self.transform == 'id': + pixdim = [1.0, 1.0, 1.0] + voxToDisplayMat = np.eye(4) + elif self.transform == 'pixdim': + pixdim = self.image.pixdim + voxToDisplayMat = np.diag([pixdim[0], pixdim[1], pixdim[2], 1.0]) + elif self.transform == 'affine': + voxToDisplayMat = self.image.voxToWorldMat + + voxToDisplayMat = np.array(voxToDisplayMat, dtype=np.float32) + displayToVoxMat = transform.invert(voxToDisplayMat) + + self.voxToDisplayMat = voxToDisplayMat.transpose() + self.displayToVoxMat = displayToVoxMat.transpose() + + # for pixdim/identity transformations, we want the world + # location (0, 0, 0) to map to voxel location (0, 0, 0) + for i in range(3): + self.voxToDisplayMat[3, i] = pixdim[i] * 0.5 + self.displayToVoxMat[3, i] = -0.5 + # When transform is changed to 'affine', enable interpolation # and, when changed to 'pixdim' or 'id', disable interpolation - def changeInterp(*a): - if self.transform == 'affine': self.interpolation = 'spline' - else: self.interpolation = 'none' - - self.addListener('transform', - 'ImageDisplay_{}'.format(id(self)), - changeInterp) + if self.transform == 'affine': self.interpolation = 'spline' + else: self.interpolation = 'none' class DisplayContext(props.HasProperties): @@ -307,6 +356,14 @@ class DisplayContext(props.HasProperties): 3D location (xyz) in the image list space. """ + + bounds = props.Bounds(ndims=3) + """This property contains the min/max values of a bounding box (in display + coordinates) which is big enough to contain all of the images in the + :attr:`images` list. This property shouid be read-only, but I don't have a + way to enforce it (yet). + """ + volume = props.Int(minval=0, maxval=0, default=0, clamped=True) """The volume property contains the currently selected volume @@ -325,12 +382,14 @@ class DisplayContext(props.HasProperties): self._imageList = imageList self._name = '{}_{}'.format(self.__class__.__name__, id(self)) + # Ensure that an ImageDisplay object exists for + # every image, and that the display bounds + # property is initialised self._imageListChanged() - self._imageListBoundsChanged() # initialise the location to be # the centre of the image world - b = imageList.bounds + b = self.bounds self.location.xyz = [ b.xlo + b.xlen / 2.0, b.ylo + b.ylen / 2.0, @@ -339,9 +398,9 @@ class DisplayContext(props.HasProperties): imageList.addListener('images', self._name, self._imageListChanged) - imageList.addListener('bounds', + self.addListener( 'bounds', self._name, - self._imageListBoundsChanged) + self._boundsChanged) self.addListener( 'volume', self._name, self._volumeChanged) @@ -361,9 +420,29 @@ class DisplayContext(props.HasProperties): # Ensure that an ImageDisplay # object exists for every image for image in self._imageList: - try: image.getAttribute('display') - except KeyError: image.setAttribute('display', ImageDisplay(image)) - + try: + display = image.getAttribute('display') + except KeyError: + display = ImageDisplay(image) + image.setAttribute('display', display) + + # Register a listener with the transform property + # of every image display so that when they change, + # we can update the display bounds. + # + # This may be called multiple times on each image, + # but it doesn't matter, as any listener which has + # previously been registered with an image will + # just be replaced by the new one here. + display.addListener( + 'transform', + self.__class__.__name__, + self._updateBounds, + overwrite=True) + + # Ensure that the bounds property is accurate + self._updateBounds() + # Limit the selectedImage property # so it cannot take a value greater # than len(imageList)-1 @@ -391,6 +470,36 @@ class DisplayContext(props.HasProperties): self.setConstraint('volume', 'maxval', 0) + def _updateBounds(self, *a): + """Called when the image list changes, or when any image display + transform is changed. Updates the :attr:`bounds` property. + """ + + if len(self._imageList) == 0: + minBounds = [0.0, 0.0, 0.0] + maxBounds = [0.0, 0.0, 0.0] + + else: + minBounds = 3 * [ sys.float_info.max] + maxBounds = 3 * [-sys.float_info.max] + + for img in self.images: + + display = img.getAttribute('display') + xform = display.voxToDisplayMat + + for ax in range(3): + + lo, hi = transform.axisBounds(img.shape, xform, ax) + + if lo < minBounds[ax]: minBounds[ax] = lo + if hi > maxBounds[ax]: maxBounds[ax] = hi + + self.bounds[:] = [minBounds[0], maxBounds[0], + minBounds[1], maxBounds[1], + minBounds[2], maxBounds[2]] + + def _volumeChanged(self, *a): """Called when the :attr:`volume` property changes. @@ -406,14 +515,12 @@ class DisplayContext(props.HasProperties): image.getAttribute('display').volume = self.volume - def _imageListBoundsChanged(self, *a): - """Called when the :attr:`fsl.data.image.ImageList.bounds` property - changes. + def _boundsChanged(self, *a): + """Called when the :attr:`bounds` property changes. Updates the constraints on the :attr:`location` property. """ - bounds = self._imageList.bounds - self.location.setLimits(0, bounds.xlo, bounds.xhi) - self.location.setLimits(1, bounds.ylo, bounds.yhi) - self.location.setLimits(2, bounds.zlo, bounds.zhi) + self.location.setLimits(0, self.bounds.xlo, self.bounds.xhi) + self.location.setLimits(1, self.bounds.ylo, self.bounds.yhi) + self.location.setLimits(2, self.bounds.zlo, self.bounds.zhi) diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py new file mode 100644 index 000000000..b564ce0e2 --- /dev/null +++ b/fsl/utils/transform.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python +# +# transform.py - Functions for applying affine transformations to +# coordinates. +# +# Author: Paul McCarthy <pauldmccarthy@gmail.com> +# +"""This module provides functions related to 3D image transformations and +spaces. +""" + +import numpy as np +import numpy.linalg as linalg +import collections + + +def invert(x): + """Inverts the given matrix. """ + return linalg.inv(x) + + +def axisBounds(shape, x, axis): + """Returns the (lo, hi) bounds of the specified axis.""" + x, y, z = shape + + x -= 0.5 + y -= 0.5 + z -= 0.5 + + points = np.zeros((8, 3), dtype=np.float32) + + 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] + + tx = transform(points, x) + + lo = tx[:, axis].min() + hi = tx[:, axis].max() + + return (lo, hi) + + +def axisLength(shape, x, axis): + """Return the length, in real world units, of the specified axis. + """ + + points = np.zeros((2, 3), dtype=np.float32) + points[:] = [-0.5, -0.5, -0.5] + points[1, axis] = shape[axis] - 0.5 + + tx = transform(points, x) + + # euclidean distance between each boundary point + return sum((tx[0, :] - tx[1, :]) ** 2) ** 0.5 + + +def transform(p, x, axes): + """Transforms the given set of points ``p`` according to the given affine + transformation ``x``. The transformed points are returned as a + :class:``numpy.float64`` array. + """ + + p = _fillPoints(p, axes) + t = np.zeros((len(p), 3), dtype=np.float64) + + x = p[:, 0] + y = p[:, 1] + z = p[:, 2] + + t[:, 0] = x * x[0, 0] + y * x[1, 0] + z * x[2, 0] + x[3, 0] + t[:, 1] = x * x[0, 1] + y * x[1, 1] + z * x[2, 1] + x[3, 1] + t[:, 2] = x * x[0, 2] + y * x[1, 2] + z * x[2, 2] + x[3, 2] + + if axes is None: axes = [0, 1, 2] + + tx = np.array(t[:, axes], dtype=np.float64) + + if tx.size == 1: return tx[0] + else: return tx + + +def _fillPoints(p, axes): + """Used by the :func:`transform` function. Turns the given array p into + a N*3 array of x,y,z coordinates. The array p may be a 1D array, or an + N*2 or N*3 array. + """ + + if not isinstance(p, collections.Iterable): p = [p] + + p = np.array(p) + + if axes is None: return p + + if not isinstance(axes, collections.Iterable): axes = [axes] + + if p.ndim == 1: + p = p.reshape((len(p), 1)) + + if p.ndim != 2: + raise ValueError('Points array must be either one or two ' + 'dimensions') + + if len(axes) != p.shape[1]: + raise ValueError('Points array shape does not match specified ' + 'number of axes') + + newp = np.zeros((len(p), 3), dtype=p.dtype) + + for i, ax in enumerate(axes): + newp[:, ax] = p[:, i] + + return newp -- GitLab