Merge branch 'master' into 'v1.0'

Version 1.0.5

See merge request !33

fsl/data/gifti.py

@@ -49,7 +49,7 @@ class GiftiSurface(mesh.TriangleMesh):
     """
 
 
-    def __init__(self, infile):
+    def __init__(self, infile, fixWinding=False):
         """Load the given GIFTI file using ``nibabel``, and extracts surface
         data using the  :func:`loadGiftiSurface` function.
 
@@ -61,7 +61,7 @@ class GiftiSurface(mesh.TriangleMesh):
 
         surfimg, vertices, indices = loadGiftiSurface(infile)
 
-        mesh.TriangleMesh.__init__(self, vertices, indices)
+        mesh.TriangleMesh.__init__(self, vertices, indices, fixWinding)
 
         name   = fslpath.removeExt(op.basename(infile), ALLOWED_EXTENSIONS)
         infile = op.abspath(infile)

fsl/data/mesh.py

@@ -27,6 +27,8 @@ import numpy   as np
 
 import six
 
+import fsl.utils.transform as transform
+
 from . import image as fslimage
 
 
@@ -34,9 +36,9 @@ log = logging.getLogger(__name__)
 
 
 class TriangleMesh(object):
-    """The ``TriangleMesh`` class represents a 3D model. A mesh is defined by
-    a collection of vertices and indices.  The indices index into the list of
-    vertices, and define a set of triangles which make the model.
+    """The ``TriangleMesh`` class represents a 3D model. A mesh is defined by a
+    collection of ``N`` vertices, and ``M`` triangles.  The triangles are
+    defined by ``(M, 3)) indices into the list of vertices.
 
 
     A ``TriangleMesh`` instance has the following attributes:
@@ -53,6 +55,12 @@ class TriangleMesh(object):
 
     ``indices``    A :meth:`M\times 3` ``numpy`` array containing
                    the vertex indices for :math:`M` triangles
+
+    ``normals``    A :math:`M\times 3` ``numpy`` array containing
+                   face normals.
+
+    ``vnormals``   A :math:`N\times 3` ``numpy`` array containing
+                   vertex normals.
     ============== ====================================================
 
 
@@ -68,16 +76,20 @@ class TriangleMesh(object):
     """
 
 
-    def __init__(self, data, indices=None):
+    def __init__(self, data, indices=None, fixWinding=False):
         """Create a ``TriangleMesh`` instance.
 
-        :arg data:    Can either be a file name, or a :math:`N\\times 3`
-                      ``numpy`` array containing vertex data. If ``data`` is
-                      a file name, it is passed to the
-                      :func:`loadVTKPolydataFile` function.
+        :arg data:       Can either be a file name, or a :math:`N\\times 3`
+                         ``numpy`` array containing vertex data. If ``data``
+                         is a file name, it is passed to the
+                         :func:`loadVTKPolydataFile` function.
+
+        :arg indices:    A list of indices into the vertex data, defining
+                         the triangles.
 
-        :arg indices: A list of indices into the vertex data, defining
-                      the triangles.
+        :arg fixWinding: Defaults to ``False``. If ``True``, the vertex
+                         winding order of every triangle is is fixed so they
+                         all have outward-facing normal vectors.
         """
 
         if isinstance(data, six.string_types):
@@ -97,12 +109,17 @@ class TriangleMesh(object):
         if indices is None:
             indices = np.arange(data.shape[0])
 
-        self.vertices     = np.array(data)
-        self.indices      = np.array(indices).reshape((-1, 3))
+        self.__vertices     = np.array(data)
+        self.__indices      = np.array(indices).reshape((-1, 3))
 
         self.__vertexData = {}
-        self.__loBounds   = self.vertices.min(axis=0)
-        self.__hiBounds   = self.vertices.max(axis=0)
+        self.__faceNormals = None
+        self.__vertNormals = None
+        self.__loBounds    = self.vertices.min(axis=0)
+        self.__hiBounds    = self.vertices.max(axis=0)
+
+        if fixWinding:
+            self.__fixWindingOrder()
 
 
     def __repr__(self):
@@ -118,6 +135,101 @@ class TriangleMesh(object):
         return self.__repr__()
 
 
+    @property
+    def vertices(self):
+        """The ``(N, 3)`` vertices of this mesh. """
+        return self.__vertices
+
+
+    @property
+    def indices(self):
+        """The ``(M, 3)`` triangles of this mesh. """
+        return self.__indices
+
+
+    def __fixWindingOrder(self):
+        """Called by :meth:`__init__` if ``fixWinding is True``.  Fixes the
+        mesh triangle winding order so that all face normals are facing
+        outwards from the centre of the mesh.
+        """
+
+        # Define a viewpoint which is
+        # far away from the mesh.
+        fnormals = self.normals
+        camera   = self.__loBounds - (self.__hiBounds - self.__loBounds)
+
+        # Find the nearest vertex
+        # to the viewpoint
+        dists = np.sqrt(np.sum((self.vertices - camera) ** 2, axis=1))
+        ivert = np.argmin(dists)
+        vert  = self.vertices[ivert]
+
+        # Pick a triangle that
+        # this vertex in and
+        # ges its face normal
+        itri = np.where(self.indices == ivert)[0][0]
+        n    = fnormals[itri, :]
+
+        # Make sure the angle between the
+        # normal, and a vector from the
+        # vertex to the camera is positive
+        # If it isn't, flip the triangle
+        # winding order.
+        if np.dot(n, transform.normalise(camera - vert)) < 0:
+            self.indices[:, [1, 2]] = self.indices[:, [2, 1]]
+            self.__faceNormals     *= -1
+
+
+    @property
+    def normals(self):
+        """A ``(M, 3)`` array containing surface normals for every
+        triangle in the mesh, normalised to unit length.
+        """
+
+        if self.__faceNormals is not None:
+            return self.__faceNormals
+
+        v0 = self.vertices[self.indices[:, 0]]
+        v1 = self.vertices[self.indices[:, 1]]
+        v2 = self.vertices[self.indices[:, 2]]
+
+        n = np.cross((v1 - v0), (v2 - v0))
+
+        self.__faceNormals = transform.normalise(n)
+
+        return self.__faceNormals
+
+
+    @property
+    def vnormals(self):
+        """A ``(N, 3)`` array containing normals for every vertex
+        in the mesh.
+        """
+        if self.__vertNormals is not None:
+            return self.__vertNormals
+
+        # per-face normals
+        fnormals = self.normals
+        vnormals = np.zeros((self.vertices.shape[0], 3), dtype=np.float)
+
+        # TODO make fast. I can't figure
+        # out how to use np.add.at to
+        # accumulate the face normals for
+        # each vertex.
+        for i in range(self.indices.shape[0]):
+
+            v0, v1, v2 = self.indices[i]
+
+            vnormals[v0, :] += fnormals[i]
+            vnormals[v1, :] += fnormals[i]
+            vnormals[v2, :] += fnormals[i]
+
+        # normalise to unit length
+        self.__vertNormals = transform.normalise(vnormals)
+
+        return self.__vertNormals
+
+
     def getBounds(self):
         """Returns a tuple of values which define a minimal bounding box that
         will contain all vertices in this ``TriangleMesh`` instance. The

fsl/utils/memoize.py

@@ -10,6 +10,8 @@ a function:
  .. autosummary::
     :nosignatures:
 
+    memoize
+    Memoize
     Instanceify
     memoizeMD5
     skipUnchanged
@@ -17,7 +19,6 @@ a function:
 
 
 import logging
-
 import hashlib
 import functools
 import six
@@ -25,20 +26,91 @@ import six
 log = logging.getLogger(__name__)
 
 
-# TODO Make this a class, and add
-#      a "clearCache" method to it.
-def memoize(func):
+def memoize(func=None):
     """Memoize the given function by the value of the input arguments.
 
+    This function simply returns a :class:`Memoize` instance.
+    """
+
+    return Memoize(func)
+
+
+class Memoize(object):
+    """Decorator which can be used to memoize a function or method. Use like
+    so::
+
+        @memoize
+        def myfunc(*a, **kwa):
+            ...
+
+        @memoize()
+        def otherfunc(*a, **kwax):
+            ...
+
+    A ``Memoize`` instance maintains a cache which contains ``{args : value}``
+    mappings, where ``args`` are the input arguments to the function, and
+    ``value`` is the value that the function returned for those arguments.
+    When a memoized function is called with arguments that are present in the
+    cache, the cached values are returned, and the function itself is not
+    called.
+
+
+    The :meth:`invalidate` method may be used to clear the internal cache.
+
+
     Note that the arguments used for memoization must be hashable, as they are
     used as keys in a dictionary.
     """
 
-    cache      = {}
-    defaultKey = '_memoize_noargs_'
 
-    def wrapper(*a, **kwa):
+    def __init__(self, *args, **kwargs):
+        """Create a ``Memoize`` object.
+        """
+
+        self.__cache      = {}
+        self.__func       = None
+        self.__defaultKey = '_memoize_noargs_'
+
+        self.__setFunction(*args, **kwargs)
+
+
+    def invalidate(self, *args, **kwargs):
+        """Clears the internal cache. If no arguments are given, the entire
+        cache is cleared. Otherwise, only the cached value for the provided
+        arguments is cleared.
+        """
+
+        if len(args) + len(kwargs) == 0:
+            self.__cache = {}
+
+        else:
+            key = self.__makeKey(*args, **kwargs)
+
+            try:
+                self.__cache.pop(key)
+            except KeyError:
+                pass
+
+
+    def __setFunction(self, *args, **kwargs):
+        """Used internally to set the memoized function. """
+
+        if self.__func is not None:
+            return False
+
+        # A no-brackets style
+        # decorator was used
+        isfunc = (len(kwargs) == 0 and len(args) == 1 and callable(args[0]))
+
+        if isfunc:
+            self.__func = args[0]
+
+        return isfunc
 
+
+    def __makeKey(self, *a, **kwa):
+        """Constructs a key for use with the cache from the given arguments.
+        """
         key = []
 
         if a   is not None: key += list(a)
@@ -48,24 +120,35 @@ def memoize(func):
         # any arguments specified - use the
         # default cache key.
         if len(key) == 0:
-            key = [defaultKey]
+            key = [self.__defaultKey]
+
+        return tuple(key)
 
-        key = tuple(key)
+
+    def __call__(self, *a, **kwa):
+        """Checks the cache against the given arguments. If a cached value
+        is present, it is returned. Otherwise the memoized function is called,
+        and its value is cached and returned.
+        """
+
+        if self.__setFunction(*a, **kwa):
+            return self
+
+        key = self.__makeKey(*a, **kwa)
 
         try:
-            result = cache[key]
+            result = self.__cache[key]
 
             log.debug(u'Retrieved from cache[{}]: {}'.format(key, result))
 
         except KeyError:
 
-            result     = func(*a, **kwa)
-            cache[key] = result
+            result            = self.__func(*a, **kwa)
+            self.__cache[key] = result
 
             log.debug(u'Adding to cache[{}]: {}'.format(key, result))
 
         return result
-    return wrapper
 
 
 def memoizeMD5(func):
@@ -143,8 +226,13 @@ def skipUnchanged(func):
             newIsArray = isinstance(value,  np.ndarray)
             isarray    = oldIsArray or newIsArray
 
-            if isarray: nochange = np.all(oldVal == value)
-            else:       nochange =        oldVal == value
+            if isarray:
+                a = np.array(oldVal, copy=False)
+                b = np.array(value,  copy=False)
+
+                nochange = (a.shape == b.shape) and np.allclose(a, b)
+            else:
+                nochange = oldVal == value
 
             if nochange:
                 return False

fsl/utils/transform.py

@@ -11,16 +11,26 @@ spaces. The following functions are provided:
    :nosignatures:
 
    transform
+   transformNormal
    scaleOffsetXform
    invert
    concat
    compose
    decompose
+   rotMatToAffine
    rotMatToAxisAngles
    axisAnglesToRotMat
    axisBounds
    flirtMatrixToSform
    sformToFlirtMatrix
+
+And a few more functions are provided for working with vectors:
+
+.. autosummary::
+   :nosignatures:
+
+   veclength
+   normalise
 """
 
 import numpy        as np
@@ -44,6 +54,29 @@ def concat(*xforms):
     return result
 
 
+def veclength(vec):
+    """Returns the length of the given vector(s).
+
+    Multiple vectors may be passed in, with a shape of ``(n, 3)``.
+    """
+    vec = np.array(vec, copy=False).reshape(-1, 3)
+    return np.sqrt(np.einsum('ij,ij->i', vec, vec))
+
+
+def normalise(vec):
+    """Normalises the given vector(s) to unit length.
+
+    Multiple vectors may be passed in, with a shape of ``(n, 3)``.
+    """
+    vec = np.array(vec, copy=False).reshape(-1, 3)
+    n   = (vec.T / veclength(vec)).T
+
+    if n.size == 3:
+        n = n[0]
+
+    return n
+
+
 def scaleOffsetXform(scales, offsets):
     """Creates and returns an affine transformation matrix which encodes
     the specified scale(s) and offset(s).
@@ -60,10 +93,12 @@ def scaleOffsetXform(scales, offsets):
     :returns:     A ``numpy.float32`` array of size :math:`4 \\times 4`.
     """
 
-    if not isinstance(scales,  collections.Sequence): scales  = [scales]
-    if not isinstance(offsets, collections.Sequence): offsets = [offsets]
-    if not isinstance(scales,  list):                 scales  = list(scales)
-    if not isinstance(offsets, list):                 offsets = list(offsets)
+    oktypes = (collections.Sequence, np.ndarray)
+
+    if not isinstance(scales,  oktypes): scales  = [scales]
+    if not isinstance(offsets, oktypes): offsets = [offsets]
+    if not isinstance(scales,  list):    scales  = list(scales)
+    if not isinstance(offsets, list):    offsets = list(offsets)
 
     lens = len(scales)
     leno = len(offsets)
@@ -131,7 +166,7 @@ def compose(scales, offsets, rotations, origin=None):
     return concat(offset, postRotate, rotate, preRotate, scale)
 
 
-def decompose(xform):
+def decompose(xform, angles=True):
     """Decomposes the given transformation matrix into separate offsets,
     scales, and rotations, according to the algorithm described in:
 
@@ -142,12 +177,17 @@ def decompose(xform):
     It is assumed that the given transform has no perspective components. Any
     shears in the affine are discarded.
 
-    :arg xform: A ``(4, 4)`` affine transformation matrix.
+    :arg xform:  A ``(4, 4)`` affine transformation matrix.
+
+    :arg angles: If ``True`` (the default), the rotations are returned
+                 as axis-angles, in radians. Otherwise, the rotation matrix
+                 is returned.
 
     :returns: The following:
                 - A sequence of three scales
                 - A sequence of three translations
-                - A sequence of three rotations, in radians
+                - A sequence of three rotations, in radians. Or, if
+                  ``angles is False``, a rotation matrix.
     """
 
     # The inline comments in the code below are taken verbatim from
@@ -216,9 +256,17 @@ def decompose(xform):
     # Finally, we need to decompose the rotation matrix into a sequence
     # of rotations about the x, y, and z axes. [This is done in the
     # rotMatToAxisAngles function]
-    rx, ry, rz = rotMatToAxisAngles(R.T)
+    if angles: rotations = rotMatToAxisAngles(R.T)
+    else:      rotations = R.T
+
+    return [sx, sy, sz], translations, rotations
+
 
-    return [sx, sy, sz], translations, [rx, ry, rz]
+def rotMatToAffine(rotmat, origin=None):
+    """Convenience function which encodes the given ``(3, 3)`` rotation
+    matrix into a ``(4, 4)`` affine.
+    """
+    return compose([1, 1, 1], [0, 0, 0], rotmat, origin)
 
 
 def rotMatToAxisAngles(rotmat):
@@ -396,24 +444,29 @@ def axisBounds(shape,
     else:      return (lo,    hi)
 
 
-def transform(p, xform, axes=None):
+def transform(p, xform, axes=None, vector=False):
     """Transforms the given set of points ``p`` according to the given affine
     transformation ``xform``.
 
 
-    :arg p:     A sequence or array of points of shape :math:`N \\times  3`.
+    :arg p:      A sequence or array of points of shape :math:`N \\times  3`.
+
+    :arg xform:  A ``(4, 4)`` affine transformation matrix with which to
+                 transform the points in ``p``.
 
-    :arg xform: An affine transformation matrix with which to transform the
-                points in ``p``.
+    :arg axes:   If you are only interested in one or two axes, and the source
+                 axes are orthogonal to the target axes (see the note below),
+                 you may pass in a 1D, ``N*1``, or ``N*2`` array as ``p``, and
+                 use this argument to specify which axis/axes that the data in