Routine to fix winding order/normal vectors on triangle meshes. I am not

convinced that this is foolproof, but it seems to work.

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:
@@ -74,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):
@@ -103,8 +109,8 @@ 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.__faceNormals = None
@@ -112,6 +118,9 @@ class TriangleMesh(object):
         self.__loBounds    = self.vertices.min(axis=0)
         self.__hiBounds    = self.vertices.max(axis=0)
 
+        if fixWinding:
+            self.__fixWindingOrder()
+
 
     def __repr__(self):
         """Returns a string representation of this ``TriangleMesh`` instance.
@@ -127,35 +136,73 @@ class TriangleMesh(object):
 
 
     @property
-    def normals(self):
-        """Returns a ``(M, 3)`` array containing normals for every triangle
-        in the mesh.
+    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.
         """
-        if self.__faceNormals is not None:
-            return self.__faceNormals
 
-        v1 = self.vertices[self.indices[:, 0]]
-        v2 = self.vertices[self.indices[:, 1]]
+        # 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
 
-        fnormals = np.cross(v1, v2)
 
-        # TODO make fast, and make sure that this actually works.
-        for i in range(self.indices.shape[0]):
+    @property
+    def normals(self):
+        """A ``(M, 3)`` array containing surface normals for every
+        triangle in the mesh, normalised to unit length.
+        """
 
-            p0, p1, p2 = self.vertices[self.indices[i], :]
-            n          = fnormals[i, :]
+        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]]
 
-            if np.dot(n, np.cross(p1 - p0, p2 - p0)) > 0:
-                fnormals[i, :] *= -1
+        n = np.cross((v1 - v0), (v2 - v0))
 
-        self.__faceNormals = fnormals
+        self.__faceNormals = transform.normalise(n)
 
         return self.__faceNormals
 
 
     @property
     def vnormals(self):
-        """Returns a ``(N, 3)`` array containing normals for every vertex
+        """A ``(N, 3)`` array containing normals for every vertex
         in the mesh.
         """
         if self.__vertNormals is not None:
@@ -163,10 +210,12 @@ class TriangleMesh(object):
 
         # per-face normals
         fnormals = self.normals
-
         vnormals = np.zeros((self.vertices.shape[0], 3), dtype=np.float)
 
-        # TODO make fast
+        # 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]
@@ -176,10 +225,7 @@ class TriangleMesh(object):
             vnormals[v2, :] += fnormals[i]
 
         # normalise to unit length
-        lens = np.sqrt(np.sum(vnormals * vnormals, axis=1))
-        vnormals = (vnormals.T / lens).T
-
-        self.__vertNormals = vnormals
+        self.__vertNormals = transform.normalise(vnormals)
 
         return self.__vertNormals