Skip to content
Snippets Groups Projects
Commit 6ca03276 authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

TriangleMesh is no more. Things are broken.

parent 1e360ab9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 294 additions and 256 deletions
fsl/data/mesh.py
+ 294
− 256
View file @ 6ca03276
...@@ -4,30 +4,40 @@ ...@@ -4,30 +4,40 @@
# #
# Author: Paul McCarthy <pauldmccarthy@gmail.com> # Author: Paul McCarthy <pauldmccarthy@gmail.com>
# #
"""This module provides the :class:`TriangleMesh` class, which represents a """This module provides the :class:`Mesh` class, which represents a
3D model made of triangles. 3D model made of triangles.
.. note:: I/O support is very limited - currently, the only supported file See also the following modules:
type is the VTK legacy file format, containing the ``POLYDATA``
dataset. the :class:`TriangleMesh` class assumes that every polygon
defined in an input file is a triangle (i.e. refers to three
vertices).
See http://www.vtk.org/wp-content/uploads/2015/04/file-formats.pdf .. autosummary::
for an overview of the VTK legacy file format.
In the future, I may or may not add support for more complex meshes. fsl.data.vtk
fsl.data.gifti
fsl.data.freesurfer
A handful of standalone functions are provided in this module, for doing various
things with meshes:
.. autosummary::
:nosignatures:
calcFaceNormals
calcVertexNormals
needsFixing
""" """
import logging import logging
import collections
import six
import deprecation
import os.path as op import os.path as op
import numpy as np import numpy as np
import six
import fsl.utils.meta as meta import fsl.utils.meta as meta
import fsl.utils.notifier as notifier
import fsl.utils.memoize as memoize import fsl.utils.memoize as memoize
import fsl.utils.transform as transform import fsl.utils.transform as transform
...@@ -37,13 +47,13 @@ from . import image as fslimage ...@@ -37,13 +47,13 @@ from . import image as fslimage
log = logging.getLogger(__name__) log = logging.getLogger(__name__)
class TriangleMesh(meta.Meta): class Mesh(notifier.Notifier, meta.Meta):
"""The ``TriangleMesh`` class represents a 3D model. A mesh is defined by a """The ``Mesh`` class represents a 3D model. A mesh is defined by a
collection of ``N`` vertices, and ``M`` triangles. The triangles are collection of ``N`` vertices, and ``M`` triangles. The triangles are
defined by ``(M, 3)`` indices into the list of vertices. defined by ``(M, 3)`` indices into the list of vertices.
A ``TriangleMesh`` instance has the following attributes: A ``Mesh`` instance has the following attributes:
============== ==================================================== ============== ====================================================
...@@ -52,134 +62,164 @@ class TriangleMesh(meta.Meta): ...@@ -52,134 +62,164 @@ class TriangleMesh(meta.Meta):
``dataSource`` Full path to the mesh file (or ``None`` if there is ``dataSource`` Full path to the mesh file (or ``None`` if there is
no file associated with this mesh). no file associated with this mesh).
``vertices`` A :math:`N\times 3` ``numpy`` array containing ``vertices`` A ``(n, 3)`` array containing the currently selected
the vertices. vertices. You can assign a vertex set key to this
attribute to change the selected vertex set.
``bounds`` The lower and upper bounds
``indices`` A :meth:`M\times 3` ``numpy`` array containing ``indices`` A ``(m, 3)`` array containing the vertex indices
the vertex indices for :math:`M` triangles for ``m`` triangles
``normals`` A :math:`M\times 3` ``numpy`` array containing ``normals`` A ``(m, 3)`` array containing face normals for the
face normals. triangles
``vnormals`` A :math:`N\times 3` ``numpy`` array containing ``vnormals`` A ``(n, 3)`` array containing vertex normals for the
vertex normals. the current vertices.
============== ==================================================== ============== ====================================================
And the following methods: **Vertex sets**
A ``Mesh`` object can be associated with multiple sets of vertices, but
only one set of triangles. Vertices can be added via the
:meth:`addVertices` method. Each vertex set must be associated with a
unique key - you can then select the current vertex set via the
:meth:`vertices` property. Most ``Mesh`` methods will raise a ``KeyError``
if you have not added any vertex sets, or selected a vertex set.
**Vertex data**
A ``Mesh`` object can store vertex-wise data. The following methods can be
used for adding/retrieving vertex data:
.. autosummary:: .. autosummary::
:nosignatures: :nosignatures:
getBounds addVertexData
loadVertexData
getVertexData getVertexData
clearVertexData clearVertexData
"""
def __init__(self, data, indices=None, fixWinding=False): **Notification**
"""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 indices: A list of indices into the vertex data, defining The ``Mesh`` class inherits from the :class:`Notifier` class. Whenever the
the triangles. ``Mesh`` vertex set is changed, a notification is emitted via the
``Notifier`` interface, with a topic of ``'vertices'``. When this occurs,
the :meth:`vertices`, :meth:`bounds`, :meth:`normals` and :attr:`vnormals`
properties will all change so that they return data specific to the newly
selected vertex set.
: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): **Metadata*
infile = data
data, lengths, indices = loadVTKPolydataFile(infile)
if np.any(lengths != 3):
raise RuntimeError('All polygons in VTK file must be '
'triangles ({})'.format(infile))
self.name = op.basename(infile) The ``Mesh`` class also inherits from the :class:`Meta` class, so
self.dataSource = infile any metadata associated with the ``Mesh`` may be added via those methods.
else:
self.name = 'TriangleMesh'
self.dataSource = None
if indices is None:
indices = np.arange(data.shape[0])
self.__vertices = np.array(data) **Geometric queries**
self.__indices = np.array(indices).reshape((-1, 3))
self.__vertexData = {}
self.__faceNormals = None
self.__vertNormals = None
self.__loBounds = self.vertices.min(axis=0)
self.__hiBounds = self.vertices.max(axis=0)
if fixWinding: If the ``trimesh`` library is present, the following methods may be used
self.__fixWindingOrder() to perform geometric queries on a mesh:
.. autosummary::
:nosignatures:
def __repr__(self): rayIntersection
"""Returns a string representation of this ``TriangleMesh`` instance. planeIntersection
""" nearestVertex
return '{}({}, {})'.format(type(self).__name__, """
self.name,
self.dataSource)
def __str__(self):
"""Returns a string representation of this ``TriangleMesh`` instance.
"""
return self.__repr__()
def __init__(self, indices, name='mesh', dataSource=None):
"""Create a ``Mesh`` instance.
Before a ``Mesh`` can be used, some vertices must be added via the
:meth:`addVertices` method.
def __fixWindingOrder(self): :arg indices: A list of indices into the vertex data, defining the
"""Called by :meth:`__init__` if ``fixWinding is True``. Fixes the mesh triangles.
mesh triangle winding order so that all face normals are facing
outwards from the centre of the mesh. :arg name: A name for this ``Mesh``.
:arg dataSource: the data source for this ``Mesh``.
""" """
# Define a viewpoint which is self.__name = name
# far away from the mesh. self.__dataSource = dataSource
fnormals = self.normals self.__indices = np.asarray(indices).reshape((-1, 3))
camera = self.__loBounds - (self.__hiBounds - self.__loBounds)
# This attribute is used to store
# Find the nearest vertex # the currently selected vertex set,
# to the viewpoint # used as a kety into all of the
dists = np.sqrt(np.sum((self.vertices - camera) ** 2, axis=1)) # dictionaries below.
ivert = np.argmin(dists) self.__selected = None
vert = self.vertices[ivert]
# Flag used to keep track of whether
# Pick a triangle that # the triangle winding order has been
# this vertex in and # "fixed" - see the addVertices method.
# ges its face normal self.__fixed = False
itri = np.where(self.indices == ivert)[0][0]
n = fnormals[itri, :] # All of these are populated
# in the addVertices method
# Make sure the angle between the self.__vertices = collections.OrderedDict()
# normal, and a vector from the self.__loBounds = collections.OrderedDict()
# vertex to the camera is positive self.__hiBounds = collections.OrderedDict()
# If it isn't, flip the triangle
# winding order. # These get populated on
if np.dot(n, transform.normalise(camera - vert)) < 0: # normals/vnormals accesses
self.indices[:, [1, 2]] = self.indices[:, [2, 1]] self.__faceNormals = collections.OrderedDict()
self.__faceNormals *= -1 self.__vertNormals = collections.OrderedDict()
# this gets populated in
# the addVertexData method
self.__vertexData = collections.OrderedDict()
# this gets populated
# in the trimesh method
self.__trimesh = collections.OrderedDict()
@property
def name(self):
"""Returns the name of this ``Mesh``. """
return self.__name
@property
def dataSource(self):
"""Returns the data source of this ``Mesh``. """
return self.__dataSource
@property @property
def vertices(self): def vertices(self):
"""The ``(N, 3)`` vertices of this mesh. """ """The ``(N, 3)`` vertices of this mesh. """
return self.__vertices return self.__vertices[self.__selected]
@property.setter
def vertices(self, key):
"""Select the current vertex set - a ``KeyError`` is raised
if no vertex set with the specified ``key`` has been added.
"""
# Force a key error if
# the key is invalid
self.__vertices[key]
self.__selected = key
@property @property
def indices(self): def indices(self):
"""The ``(M, 3)`` triangles of this mesh. """ """The ``(M, 3)`` triangles of this mesh. """
return self.__indices return self.__indices[self.__selected]
@property @property
...@@ -188,18 +228,16 @@ class TriangleMesh(meta.Meta): ...@@ -188,18 +228,16 @@ class TriangleMesh(meta.Meta):
triangle in the mesh, normalised to unit length. triangle in the mesh, normalised to unit length.
""" """
if self.__faceNormals is not None: selected = self.__selected
return self.__faceNormals indices = self.__indices
vertices = self.__vertices[selected]
v0 = self.vertices[self.indices[:, 0]] fnormals = self.__faceNormals.get(selected, None)
v1 = self.vertices[self.indices[:, 1]]
v2 = self.vertices[self.indices[:, 2]]
n = np.cross((v1 - v0), (v2 - v0))
self.__faceNormals = transform.normalise(n) if fnormals is None:
fnormals = calcFaceNormals(vertices, indices)
self.__faceNormals[selected] = fnormals
return self.__faceNormals return fnormals
@property @property
...@@ -207,89 +245,98 @@ class TriangleMesh(meta.Meta): ...@@ -207,89 +245,98 @@ class TriangleMesh(meta.Meta):
"""A ``(N, 3)`` array containing normals for every vertex """A ``(N, 3)`` array containing normals for every vertex
in the mesh. in the mesh.
""" """
if self.__vertNormals is not None:
return self.__vertNormals
# per-face normals selected = self.__selected
fnormals = self.normals indices = self.__indices
vnormals = np.zeros((self.vertices.shape[0], 3), dtype=np.float) vertices = self.__vertices[selected]
vnormals = self.__vertNormals.get(selected, None)
# TODO make fast. I can't figure if vnormals is None:
# out how to use np.add.at to vnormals = calcVertexNormals(vertices, indices, self.normals)
# accumulate the face normals for self.__vertNormals[selected] = vnormals
# each vertex.
for i in range(self.indices.shape[0]):
v0, v1, v2 = self.indices[i] return vnormals
vnormals[v0, :] += fnormals[i]
vnormals[v1, :] += fnormals[i]
vnormals[v2, :] += fnormals[i]
# normalise to unit length @property
self.__vertNormals = transform.normalise(vnormals) def bounds(self):
return self.__vertNormals
def getBounds(self):
"""Returns a tuple of values which define a minimal bounding box that """Returns a tuple of values which define a minimal bounding box that
will contain all vertices in this ``TriangleMesh`` instance. The will contain all of the currently selected vertices in this
bounding box is arranged like so: ``Mesh`` instance. The bounding box is arranged like so:
``((xlow, ylow, zlow), (xhigh, yhigh, zhigh))`` ``((xlow, ylow, zlow), (xhigh, yhigh, zhigh))``
""" """
return (self.__loBounds, self.__hiBounds)
lo = self.__loBounds[self.__selected]
hi = self.__hiBounds[self.__selected]
return lo, hi
def loadVertexData(self, dataSource, vertexData=None):
"""Attempts to load scalar data associated with each vertex of this
``TriangleMesh`` from the given ``dataSource``. The data is returned,
and also stored in an internal cache so it can be retrieved later
via the :meth:`getVertexData` method.
This method may be overridden by sub-classes. def addVertices(self, vertices, key, select=True, fixWinding=False):
"""Adds a set of vertices to this ``Mesh``.
:arg dataSource: Path to the vertex data to load :arg vertices: A `(n, 3)` array containing ``n`` vertices, compatible
:arg vertexData: The vertex data itself, if it has already been with the indices specified in :meth:`__init__`.
loaded.
:returns: A ``(M, N)``) array, which contains ``N`` data points :arg key: A key for this vertex set.
for ``M`` vertices.
:arg select: If ``True`` (the default), this vertex set is
made the currently selected vertex set.
: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.
""" """
nvertices = self.vertices.shape[0] vertices = np.asarray(vertices)
lo = vertices.min(axis=0)
hi = vertices.max(axis=0)
# Currently only white-space delimited self.__vertices[key] = vertices
# text files are supported self.__loBounds[key] = lo
if vertexData is None: self.__hiBounds[key] = hi
vertexData = np.loadtxt(dataSource)
vertexData.reshape(nvertices, -1)
if vertexData.shape[0] != nvertices: if select:
raise ValueError('Incompatible size: {}'.format(dataSource)) self.vertices = key
self.__vertexData[dataSource] = vertexData # indices already fixed?
if fixWinding and (not self.__fixed):
indices = self.indices
normals = self.normals
needsFix = needsFixing(vertices, indices, normals, lo, hi)
self.__fixed = True
return vertexData # See needsFixing documentation
if needsFix:
indices[:, [1, 2]] = indices[:, [2, 1]]
for k, fn in self.__faceNormals.items():
self.__faceNormals[k] = fn * -1
def addVertexData(self, key, vdata):
"""Adds a vertex-wise data set to the ``Mesh``. It can be retrieved
by passing the specified ``key`` to the :meth:`getVertexData` method.
"""
self.__vertexData[key] = vdata
def getVertexData(self, dataSource): def getVertexData(self, key):
"""Returns the vertex data for the given ``dataSource`` from the """Returns the vertex data for the given ``key`` from the
internal vertex data cache. If the given ``dataSource`` is not internal vertex data cache. If there is no vertex data iwth the
in the cache, it is loaded via :meth:`loadVertexData`. given key, a ``KeyError`` is raised.
""" """
try: return self.__vertexData[dataSource] return self.__vertexData[key]
except KeyError: return self.loadVertexData(dataSource)
def clearVertexData(self): def clearVertexData(self):
"""Clears the internal vertex data cache - see the """Clears the internal vertex data cache - see the
:meth:`loadVertexData` and :meth:`getVertexData` methods. :meth:`addVertexData` and :meth:`getVertexData` methods.
""" """
self.__vertexData = collections.OrderedDict()
self.__vertexData = {}
@memoize.Instanceify(memoize.memoize) @memoize.Instanceify(memoize.memoize)
...@@ -298,7 +345,8 @@ class TriangleMesh(meta.Meta): ...@@ -298,7 +345,8 @@ class TriangleMesh(meta.Meta):
geometric operations on the mesh. geometric operations on the mesh.
If the ``trimesh`` or ``rtree`` libraries are not available, this If the ``trimesh`` or ``rtree`` libraries are not available, this
function returns ``None`` function returns ``None``, and none of the geometric query methods
will do anything.
""" """
# trimesh is an optional dependency - rtree # trimesh is an optional dependency - rtree
...@@ -313,15 +361,17 @@ class TriangleMesh(meta.Meta): ...@@ -313,15 +361,17 @@ class TriangleMesh(meta.Meta):
log.warning('trimesh is not available') log.warning('trimesh is not available')
return None return None
if hasattr(self, '__trimesh'): tm = self.__trimesh.get(self.__selected, None)
return self.__trimesh
if tm is None:
tm = trimesh.Trimesh(self.vertices,
self.indices,
process=False,
validate=False)
self.__trimesh = trimesh.Trimesh(self.__vertices, self.__trimesh[self.__selected] = tm
self.__indices,
process=False,
validate=False)
return self.__trimesh return tm
def rayIntersection(self, origins, directions, vertices=False): def rayIntersection(self, origins, directions, vertices=False):
...@@ -453,107 +503,95 @@ class TriangleMesh(meta.Meta): ...@@ -453,107 +503,95 @@ class TriangleMesh(meta.Meta):
return lines, faces, dists return lines, faces, dists
def calcFaceNormals(vertices, indices):
"""Calculates face normals for the mesh described by ``vertices`` and
``indices``.
ALLOWED_EXTENSIONS = ['.vtk'] :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
"""A list of file extensions which could contain :class:`TriangleMesh` data. :arg indices: A ``(m, 3)`` array containing the mesh triangles.
""" :returns: A ``(m, 3)`` array containing normals for every triangle in
the mesh.
EXTENSION_DESCRIPTIONS = ['VTK polygon model file']
"""A description for each of the extensions in :data:`ALLOWED_EXTENSIONS`."""
def loadVTKPolydataFile(infile):
"""Loads a vtk legacy file containing a ``POLYDATA`` data set.
:arg infile: Name of a file to load from.
:returns: a tuple containing three values:
- A :math:`N\\times 3` ``numpy`` array containing :math:`N`
vertices.
- A 1D ``numpy`` array containing the lengths of each polygon.
- A 1D ``numpy`` array containing the vertex indices for all
polygons.
""" """
lines = None v0 = vertices[indices[:, 0]]
v1 = vertices[indices[:, 1]]
v2 = vertices[indices[:, 2]]
with open(infile, 'rt') as f: fnormals = np.cross((v1 - v0), (v2 - v0))
lines = f.readlines() fnormals = transform.normalise(fnormals)
lines = [l.strip() for l in lines] return fnormals
if lines[3] != 'DATASET POLYDATA':
raise ValueError('Only the POLYDATA data type is supported')
nVertices = int(lines[4].split()[1]) def calcVertexNormals(vertices, indices, fnormals):
nPolygons = int(lines[5 + nVertices].split()[1]) """Calculates vertex normals for the mesh described by ``vertices``
nIndices = int(lines[5 + nVertices].split()[2]) - nPolygons and ``indices``.
vertices = np.zeros((nVertices, 3), dtype=np.float32) :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
polygonLengths = np.zeros( nPolygons, dtype=np.uint32) :arg indices: A ``(m, 3)`` array containing the mesh triangles.
indices = np.zeros( nIndices, dtype=np.uint32) :arg fnormals: A ``(m, 3)`` array containing the face/triangle normals.
:returns: A ``(n, 3)`` array containing normals for every vertex in
for i in range(nVertices): the mesh.
vertLine = lines[i + 5] """
vertices[i, :] = [float(w) for w in vertLine.split()]
indexOffset = 0
for i in range(nPolygons):
polyLine = lines[6 + nVertices + i].split()
polygonLengths[i] = int(polyLine[0])
start = indexOffset vnormals = np.zeros((vertices.shape[0], 3), dtype=np.float)
end = indexOffset + polygonLengths[i]
indices[start:end] = [int(w) for w in polyLine[1:]]
indexOffset += polygonLengths[i] # 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(indices.shape[0]):
return vertices, polygonLengths, indices v0, v1, v2 = indices[i]
vnormals[v0, :] += fnormals[i]
vnormals[v1, :] += fnormals[i]
vnormals[v2, :] += fnormals[i]
def getFIRSTPrefix(modelfile): # normalise to unit length
"""If the given ``vtk`` file was generated by `FIRST return transform.normalise(vnormals)
<https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FIRST>`_, this function
will return the file prefix. Otherwise a ``ValueError`` will be
raised.
"""
if not modelfile.endswith('first.vtk'):
raise ValueError('Not a first vtk file: {}'.format(modelfile))
modelfile = op.basename(modelfile) def needsFixing(vertices, indices, fnormals, loBounds, hiBounds):
prefix = modelfile.split('-') """Determines whether the triangle winding order, for the mesh described by
prefix = '-'.join(prefix[:-1]) ``vertices`` and ``indices``, needs to be flipped.
return prefix If this function returns ``True``, the given ``indices`` and ``fnormals``
need to be adjusted so that all face normals are facing outwards from the
centre of the mesh. The necessary adjustments are as follows::
indices[:, [1, 2]] = indices[:, [2, 1]]
fnormals = fnormals * -1
def findReferenceImage(modelfile): :arg vertices: A ``(n, 3)`` array containing the mesh vertices.
"""Given a ``vtk`` file, attempts to find a corresponding ``NIFTI`` :arg indices: A ``(m, 3)`` array containing the mesh triangles.
image file. Return the path to the image, or ``None`` if no image was :arg fnormals: A ``(m, 3)`` array containing the face/triangle normals.
found. :arg loBounds: A ``(3, )`` array contaning the low vertex bounds.
:arg hiBounds: A ``(3, )`` array contaning the high vertex bounds.
Currently this function will only return an image for ``vtk`` files :returns: ``True`` if the ``indices`` and ``fnormals`` need to be
generated by `FIRST <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FIRST>`_. adjusted, ``False`` otherwise.
""" """
try: # Define a viewpoint which is
# far away from the mesh.
dirname = op.dirname(modelfile) camera = loBounds - (hiBounds - loBounds)
prefixes = [getFIRSTPrefix(modelfile)]
except ValueError: # Find the nearest vertex
return None # to the viewpoint
dists = np.sqrt(np.sum((vertices - camera) ** 2, axis=1))
if prefixes[0].endswith('_first'): ivert = np.argmin(dists)
prefixes.append(prefixes[0][:-6]) vert = vertices[ivert]
for p in prefixes: # Pick a triangle that
try: # this vertex is in and
return fslimage.addExt(op.join(dirname, p), mustExist=True) # ges its face normal
except fslimage.PathError: itri = np.where(indices == ivert)[0][0]
continue n = fnormals[itri, :]
return None # Make sure the angle between the
# normal, and a vector from the
# vertex to the camera is positive
# If it isn't, we need to flip the
# triangle winding order.
return np.dot(n, transform.normalise(camera - vert)) < 0
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment