Fix boundary_tri_index (#832)

* Fix boundary_tri_index

As far as I can tell this was actually broken because it
didn't properly account for repeated but "swapped" edges meaning
that all edges were counted as unique. From manual testing it
appears to be that now this is actually correct. It also meant
changing the ordering of the returned edges to make iterating
easier.

* Fix trimesh unit tests due to edge reordering

The ground truth values are permuted based on the change in
ordering. Also use np.testing.assert_allclose

* Add test for boundary_tri_index()

Shows that this method works as expected

menpo/shape/mesh/base.py

@@ -1,12 +1,11 @@
 # coding=utf-8
-from collections import Counter
-import numpy as np
 from warnings import warn
 
+import numpy as np
+
+from .normals import compute_vertex_normals, compute_face_normals
 from .. import PointCloud
 from ..adjacency import mask_adjacency_array, reindex_adjacency_array
-from .normals import compute_vertex_normals, compute_face_normals
-
 
 Delaunay = None  # expensive, from scipy.spatial
 
@@ -133,6 +132,7 @@ class TriMesh(PointCloud):
         Any trilist will also not be copied.
         In general this should only be used if you know what you are doing.
     """
+
     def __init__(self, points, trilist=None, copy=True):
         super(TriMesh, self).__init__(points, copy=copy)
         if trilist is None:
@@ -420,7 +420,14 @@ def mean_tri_area(self):
         return np.mean(self.tri_areas())
 
     def boundary_tri_index(self):
-        r"""Boolean index into triangles that are at the edge of the TriMesh
+        r"""Boolean index into triangles that are at the edge of the TriMesh.
+        The boundary vertices can be visualized as follows
+        ::
+
+            tri_mask = mesh.boundary_tri_index()
+            boundary_points = mesh.points[mesh.trilist[tri_mask].ravel()]
+            pc = menpo.shape.PointCloud(boundary_points)
+            pc.view()
 
         Returns
         -------
@@ -429,20 +436,35 @@ def boundary_tri_index(self):
             also an edge of another triangle (and so this triangle exists on
             the boundary of the TriMesh)
         """
-        trilist = self.trilist
-        # Get a sorted list of edge pairs
-        edge_pairs = np.sort(np.vstack((trilist[:, [0, 1]],
-                                        trilist[:, [0, 2]],
-                                        trilist[:, [1, 2]])))
-
-        # convert to a tuple per edge pair
-        edges = [tuple(x) for x in edge_pairs]
-        # count the occurrences of the ordered edge pairs - edge pairs that
-        # occur once are at the edge of the whole mesh
-        mesh_edges = (e for e, i in Counter(edges).items() if i == 1)
-        # index back into the edges to find which triangles contain these edges
-        return np.array(list(set(edges.index(e) % trilist.shape[0]
-                                 for e in mesh_edges)))
+        # Compute the edge indices so that we can find duplicated edges
+        edge_indices = self.edge_indices()
+        # Compute the triangle indices and repeat them so that when we loop
+        # over the edges we get the correct triangle index per edge
+        # (e.g. [0, 0, 0, 1, 1, 1, ...])
+        tri_indices = np.arange(self.trilist.shape[0]).repeat(3)
+
+        # Loop over the edges to find the "lonely" triangles that have an edge
+        # that isn't shared with another triangle. Due to the definition of a
+        # triangle and the careful ordering chosen above, each edge will be
+        # seen either exactly once or exactly twice.
+        # Note that some triangles may appear more than once as it's possible
+        # for a triangle to only share one edge with the rest of the mesh (so
+        # it would have two "lonely" edges
+        lonely_triangles = {}
+        for edge, t_i in zip(edge_indices, tri_indices):
+            # Sorted the edge indices since we may see an edge (0, 1) and then
+            # see it again as (1, 0) when in fact that is the same edge
+            sorted_edge = tuple(sorted(edge))
+            if sorted_edge not in lonely_triangles:
+                lonely_triangles[sorted_edge] = t_i
+            else:
+                # If we've already seen the edge the we will never see it again
+                # so we can just remove it from the candidate set
+                del lonely_triangles[sorted_edge]
+
+        mask = np.zeros(self.n_tris, dtype=np.bool)
+        mask[np.array(list(lonely_triangles.values()))] = True
+        return mask
 
     def edge_vectors(self):
         r"""A vector of edges of each triangle face.
@@ -455,14 +477,14 @@ def edge_vectors(self):
         -------
         edges : ``(n_tris * 3, n_dims)`` `ndarray`
             For each triangle (ABC), returns the edge vectors AB, BC, CA. All
-            edges are concatenated for a total of ``n_tris * 3`` edges. The
-            ordering is done so that all AB vectors are first in the returned
-            list, followed by BC, then CA.
+            edges are concatenated for a total of ``n_tris * 3`` edges.
+            The ordering is done so that each triangle is returned in order
+            e.g. [AB_1, BC_1, CA_1, AB_2, BC_2, CA_2, ...]
         """
         t = self.points[self.trilist]
-        return np.vstack((t[:, 1] - t[:, 0],
+        return np.hstack((t[:, 1] - t[:, 0],
                           t[:, 2] - t[:, 1],
-                          t[:, 2] - t[:, 0]))
+                          t[:, 2] - t[:, 0])).reshape(-1, 2)
 
     def edge_indices(self):
         r"""An unordered index into points that rebuilds the edges of this
@@ -476,14 +498,15 @@ def edge_indices(self):
         -------
         edge_indices : ``(n_tris * 3, 2)`` `ndarray`
             For each triangle (ABC), returns the pair of point indices that
-            rebuild AB, AC, BC. All edge indices are concatenated for a total
-            of ``n_tris * 3`` edge_indices. The ordering is done so that all
-            AB vectors are first in the returned list, followed by BC, then CA.
+            rebuild AB, BC, CA. All edge indices are concatenated for a total
+            of ``n_tris * 3`` edge_indices. The ordering is done so that each
+            triangle is returned in order
+            e.g. [AB_1, BC_1, CA_1, AB_2, BC_2, CA_2, ...]
         """
         tl = self.trilist
-        return np.vstack((tl[:, [0, 1]],
+        return np.hstack((tl[:, [0, 1]],
                           tl[:, [1, 2]],
-                          tl[:, [2, 0]]))
+                          tl[:, [2, 0]])).reshape(-1, 2)
 
     def unique_edge_indices(self):
         r"""An unordered index into points that rebuilds the unique edges of
@@ -719,26 +742,26 @@ def _view_2d(self, figure_id=None, new_figure=False, image_view=True,
         return PointGraphViewer2d(
             figure_id, new_figure, self.points,
             trilist_to_adjacency_array(self.trilist)).render(
-                image_view=image_view, render_lines=render_lines,
-                line_colour=line_colour, line_style=line_style,
-                line_width=line_width, render_markers=render_markers,
-                marker_style=marker_style, marker_size=marker_size,
-                marker_face_colour=marker_face_colour,
-                marker_edge_colour=marker_edge_colour,
-                marker_edge_width=marker_edge_width,
-                render_numbering=render_numbering,
-                numbers_horizontal_align=numbers_horizontal_align,
-                numbers_vertical_align=numbers_vertical_align,
-                numbers_font_name=numbers_font_name,
-                numbers_font_size=numbers_font_size,
-                numbers_font_style=numbers_font_style,
-                numbers_font_weight=numbers_font_weight,
-                numbers_font_colour=numbers_font_colour, render_axes=render_axes,
-                axes_font_name=axes_font_name, axes_font_size=axes_font_size,
-                axes_font_style=axes_font_style,
-                axes_font_weight=axes_font_weight, axes_x_limits=axes_x_limits,
-                axes_y_limits=axes_y_limits, axes_x_ticks=axes_x_ticks,
-                axes_y_ticks=axes_y_ticks, figure_size=figure_size, label=label)
+            image_view=image_view, render_lines=render_lines,
+            line_colour=line_colour, line_style=line_style,
+            line_width=line_width, render_markers=render_markers,
+            marker_style=marker_style, marker_size=marker_size,
+            marker_face_colour=marker_face_colour,
+            marker_edge_colour=marker_edge_colour,
+            marker_edge_width=marker_edge_width,
+            render_numbering=render_numbering,
+            numbers_horizontal_align=numbers_horizontal_align,
+            numbers_vertical_align=numbers_vertical_align,
+            numbers_font_name=numbers_font_name,
+            numbers_font_size=numbers_font_size,
+            numbers_font_style=numbers_font_style,
+            numbers_font_weight=numbers_font_weight,
+            numbers_font_colour=numbers_font_colour, render_axes=render_axes,
+            axes_font_name=axes_font_name, axes_font_size=axes_font_size,
+            axes_font_style=axes_font_style,
+            axes_font_weight=axes_font_weight, axes_x_limits=axes_x_limits,
+            axes_y_limits=axes_y_limits, axes_x_ticks=axes_x_ticks,
+            axes_y_ticks=axes_y_ticks, figure_size=figure_size, label=label)
 
     def _view_landmarks_2d(self, group=None, with_labels=None,
                            without_labels=None, figure_id=None,

menpo/shape/mesh/test/test_trimesh_utils.py

@@ -1,5 +1,5 @@
-from menpo.shape import TriMesh
 import numpy as np
+from menpo.shape import TriMesh
 
 
 def utils_mesh():
@@ -13,17 +13,17 @@ def utils_mesh():
 
 
 gt_edge_indices = np.array([[0, 1],
-                            [0, 2],
                             [1, 2],
-                            [2, 3],
                             [2, 0],
+                            [0, 2],
+                            [2, 3],
                             [3, 0]])
 
 gt_edge_vectors = np.array([[1, 0],
-                            [1, 1],
                             [0, 1],
-                            [-1, 0],
                             [1, 1],
+                            [1, 1],
+                            [-1, 0],
                             [0, 1]])
 
 gt_unique_edge_vectors = np.array([[1, 0],
@@ -32,36 +32,37 @@ def utils_mesh():
                                    [0, 1],
                                    [-1, 0]])
 
-gt_edge_lengths = np.array([1., 1.41421356, 1., 1., 1.41421356, 1.])
+gt_edge_lengths = np.array([1., 1., 1.41421356, 1.41421356, 1., 1.])
 
 gt_unique_edge_lengths = np.array([1., 1.41421356, 1., 1., 1.])
 
 gt_tri_areas = np.array([0.5, 0.5])
 
 
 def test_edge_indices():
-    assert np.all(utils_mesh().edge_indices() == gt_edge_indices)
+    np.testing.assert_allclose(utils_mesh().edge_indices(), gt_edge_indices)
 
 
 def test_edge_vectors():
-    assert np.all(utils_mesh().edge_vectors() == gt_edge_vectors)
+    np.testing.assert_allclose(utils_mesh().edge_vectors(), gt_edge_vectors)
 
 
 def test_unique_edge_vectors():
-    assert np.all(utils_mesh().unique_edge_vectors() == gt_unique_edge_vectors)
+    np.testing.assert_allclose(utils_mesh().unique_edge_vectors(),
+                               gt_unique_edge_vectors)
 
 
 def test_edge_lengths():
     assert np.allclose(utils_mesh().edge_lengths(), gt_edge_lengths)
 
 
 def test_unique_edge_lengths():
-    assert np.allclose(utils_mesh().unique_edge_lengths(),
-                       gt_unique_edge_lengths)
+    np.testing.assert_allclose(utils_mesh().unique_edge_lengths(),
+                               gt_unique_edge_lengths)
 
 
 def test_tri_areas():
-    assert np.allclose(utils_mesh().tri_areas(), gt_tri_areas)
+    np.testing.assert_allclose(utils_mesh().tri_areas(), gt_tri_areas)
 
 
 def test_2d_trimesh_2d_positive_areas():
@@ -72,14 +73,14 @@ def test_2d_trimesh_2d_positive_areas():
 
 
 def test_mean_tri_area():
-    assert utils_mesh().mean_tri_area() == 0.5
+    assert np.isclose(utils_mesh().mean_tri_area(), 0.5)
 
 
 def test_mean_edge_length():
-    assert np.allclose(utils_mesh().mean_edge_length(),
-                       np.mean(gt_unique_edge_lengths))
+    np.testing.assert_allclose(utils_mesh().mean_edge_length(),
+                               np.mean(gt_unique_edge_lengths))
 
 
 def test_mean_edge_length_not_unique():
-    assert np.allclose(utils_mesh().mean_edge_length(unique=False),
-                       np.mean(gt_edge_lengths))
+    np.testing.assert_allclose(utils_mesh().mean_edge_length(unique=False),
+                               np.mean(gt_edge_lengths))

menpo/shape/mesh/test/test_trimish.py

@@ -1,6 +1,8 @@
 import warnings
+
 import numpy as np
 from numpy.testing import assert_allclose
+
 from menpo.image import Image, MaskedImage
 from menpo.shape import TriMesh, TexturedTriMesh, ColouredTriMesh
 from menpo.testing import is_same_array
@@ -246,7 +248,7 @@ def test_trimesh_n_dims():
     trilist = np.array([[0, 1, 3],
                         [1, 2, 3]])
     trimesh = TriMesh(points, trilist=trilist)
-    assert(trimesh.n_dims == 3)
+    assert (trimesh.n_dims == 3)
 
 
 def test_trimesh_n_points():
@@ -257,7 +259,7 @@ def test_trimesh_n_points():
     trilist = np.array([[0, 1, 3],
                         [1, 2, 3]])
     trimesh = TriMesh(points, trilist=trilist)
-    assert(trimesh.n_points == 4)
+    assert (trimesh.n_points == 4)
 
 
 def test_trimesh_n_tris():
@@ -268,7 +270,7 @@ def test_trimesh_n_tris():
     trilist = np.array([[0, 1, 3],
                         [1, 2, 3]])
     trimesh = TriMesh(points, trilist=trilist)
-    assert(trimesh.n_tris == 2)
+    assert (trimesh.n_tris == 2)
 
 
 def test_trimesh_from_tri_mask():
@@ -281,8 +283,8 @@ def test_trimesh_from_tri_mask():
     mask = np.zeros(2, dtype=np.bool)
     mask[0] = True
     trimesh = TriMesh(points, trilist=trilist).from_tri_mask(mask)
-    assert(trimesh.n_tris == 1)
-    assert(trimesh.n_points == 3)
+    assert (trimesh.n_tris == 1)
+    assert (trimesh.n_points == 3)
     assert_allclose(trimesh.points, points[trilist[0]])
 
 
@@ -293,7 +295,7 @@ def test_trimesh_face_normals():
                        [0.0, 1.0, 0.0]])
     trilist = np.array([[0, 1, 3],
                         [1, 2, 3]])
-    expected_normals = np.array([[-np.sqrt(3)/3, -np.sqrt(3)/3, np.sqrt(3)/3],
+    expected_normals = np.array([[-np.sqrt(3) / 3, -np.sqrt(3) / 3, np.sqrt(3) / 3],
                                  [-0, -0, 1]])
     trimesh = TriMesh(points, trilist=trilist)
     face_normals = trimesh.tri_normals()
@@ -314,14 +316,35 @@ def test_trimesh_vertex_normals():
     # 0 and 2 are the corner of the triangles and so the maintain the
     # face normals. The other two are the re-normalized vertices:
     # normalize(n0 + n2)
-    expected_normals = np.array([[-np.sqrt(3)/3, -np.sqrt(3)/3, np.sqrt(3)/3],
-                                 [-0.32505758,  -0.32505758, 0.88807383],
+    expected_normals = np.array([[-np.sqrt(3) / 3, -np.sqrt(3) / 3, np.sqrt(3) / 3],
+                                 [-0.32505758, -0.32505758, 0.88807383],
                                  [0, 0, 1],
-                                 [-0.32505758,  -0.32505758, 0.88807383]])
+                                 [-0.32505758, -0.32505758, 0.88807383]])
     trimesh = TriMesh(points, trilist)
     vertex_normals = trimesh.vertex_normals()
     assert_allclose(vertex_normals, expected_normals)
 
     trimesh.points = trimesh.points.astype(np.float32)
     vertex_normals = trimesh.vertex_normals()
     assert_allclose(vertex_normals, expected_normals)
+
+
+def test_trimesh_boundary_tri_index():
+    points = np.array([[0.0, 0.0, 0.0],
+                       [0.5, 0.5, 0.0],
+                       [0.0, 0.5, 0.0],
+                       [-0.5, 0.5, 0.0],
+                       [-0.5, -0.5, 0.0],
+                       [0.5, -0.5, 0.0],
+                       [0.0, -1.0, 0.0]])
+    trilist = np.array([[0, 2, 3],
+                        [2, 0, 1],
+                        [4, 0, 3],
+                        [0, 5, 1],
+                        [4, 5, 0],
+                        [5, 4, 6]])
+    trimesh = TriMesh(points, trilist)
+    boundary_tri_index = trimesh.boundary_tri_index()
+    # The "middle" triangle is [4, 5, 0] which is surrounded on all sides
+    # [5, 4, 6] has two edges that have no neighbours
+    assert_allclose(boundary_tri_index, [True, True, True, True, False, True])