Commit cb48872f authored by Paul McCarthy's avatar Paul McCarthy 🚵
Browse files

Merge branch 'master' into 'master'

Master

See merge request !32
parents 97136a3f 2f1eb997
Pipeline #804 passed with stages
in 24 minutes and 44 seconds
......@@ -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)
......
......@@ -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
......
......@@ -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
......
......@@ -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
``p`` correspond to.
: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
``p`` correspond to.
:arg vector: Defaults to ``False``. If ``True``, the points are treated
as vectors - the translation component of the transformation
is not applied. If you set this flag, you pass in a ``(3, 3)``
transformation matrix.
:returns: The points in ``p``, transformed by ``xform``, as a ``numpy``
array with the same data type as the input.
:returns: The points in ``p``, transformed by ``xform``, as a ``numpy``
array with the same data type as the input.
.. note:: The ``axes`` argument should only be used if the source
......@@ -426,7 +479,10 @@ def transform(p, xform, axes=None):
"""
p = _fillPoints(p, axes)
t = np.dot(xform[:3, :3], p.T).T + xform[:3, 3]
t = np.dot(xform[:3, :3], p.T).T
if not vector:
t = t + xform[:3, 3]
if axes is not None:
t = t[:, axes]
......@@ -435,6 +491,14 @@ def transform(p, xform, axes=None):
else: return t
def transformNormal(p, xform, axes=None):
"""Transforms the given point(s), under the assumption that they
are normal vectors. In this case, the points are transformed by
``invert(xform[:3, :3]).T``.
"""
return transform(p, invert(xform[:3, :3]).T, axes, vector=True)
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,
......
......@@ -723,9 +723,9 @@ def _test_Image_changeData(imgtype):
assert img.saveState
assert np.all(np.isclose(img.dataRange, (dmin, dmax)))
randval = dmin + np.random.random() * drange
rx, ry, rz = randvox()
# random value within the existing data range
randval = dmin + np.random.random() * drange
rx, ry, rz = randvox()
img[rx, ry, rz] = randval
assert np.isclose(img[rx, ry, rz], randval)
......@@ -738,23 +738,29 @@ def _test_Image_changeData(imgtype):
newdmin = dmin - 100
newdmax = dmax + 100
rx, ry, rz = randvox()
img[rx, ry, rz] = newdmin
# random value below the data range
minx, miny, minz = randvox()
img[minx, miny, minz] = newdmin
assert notified.get('data', False)
assert notified.get('dataRange', False)
assert np.isclose(img[rx, ry, rz], newdmin)
assert np.isclose(img[minx, miny, minz], newdmin)
assert np.all(np.isclose(img.dataRange, (newdmin, dmax)))
notified.pop('data')
notified.pop('dataRange')
rx, ry, rz = randvox()
img[rx, ry, rz] = newdmax
# random value above the data range
maxx, maxy, maxz = randvox()
while (maxx, maxy, maxz) == (minx, miny, minz):
maxx, maxy, maxz = randvox()
img[maxx, maxy, maxz] = newdmax
assert notified.get('data', False)
assert notified.get('dataRange', False)
assert np.isclose(img[rx, ry, rz], newdmax)
assert np.isclose(img[maxx, maxy, maxz], newdmax)
assert np.all(np.isclose(img.dataRange, (newdmin, newdmax)))
finally:
......
......@@ -5,6 +5,7 @@
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
import collections
import six
import numpy as np
......@@ -50,6 +51,72 @@ def test_memoize():
assert timesCalled[0] == 7
def test_memoize_create():
timesCalled = {
'without_brackets' : 0,
'with_brackets' : 0
}
@memoize.memoize