Repository: mikedh/trimesh
Pull request #455 (2019-06-12 21:16:41)

add return_counts to unique_bincount, vertex face lookup
Makes this function more like np.unique. Plus the counts are already calculated as part of bincount.
Vertex face lookup addresses #232 using the faces_sparse array. Also falls back to a loop if sparse matrix is not available.

coveralls: 
[![Coverage Status](https://coveralls.io/builds/23984052/badge)](https://coveralls.io/builds/23984052)

Coverage decreased (-0.05%) to 82.561% when pulling **11619a811409cf768f8bdfdd7d31932abb96bc9d on BerkeleyAutomation:master** into **5aab6ade29013341868828fd6612676595f943ac on mikedh:master**.


Patch:
diff --git a/trimesh/base.py b/trimesh/base.py
index d7911cba..9d598edf 100644
--- a/trimesh/base.py
+++ b/trimesh/base.py
@@ -458,6 +458,27 @@ class Trimesh(Geometry):
             sparse=self.faces_sparse)
         return vertex_normals
 
+    @caching.cache_decorator
+    def vertex_faces(self):
+        """
+        A representation of the face indices that correspond to each vertex.
+
+        Returns
+        ----------
+        vertex_faces : (n,m) int
+          Each row contains the face indices that correspond to the given vertex,
+          padded with -1 up to the max number of faces corresponding to any one vertex
+          Where n == len(self.vertices), m == max number of faces for a single vertex
+        """
+        # make sure we have faces_sparse
+        assert hasattr(self.faces_sparse, 'dot')
+        vertex_faces = geometry.vertex_face_indices(
+            vertex_count=len(self.vertices),
+            faces=self.faces,
+            sparse=self.faces_sparse
+        )
+        return vertex_faces
+
     @vertex_normals.setter
     def vertex_normals(self, values):
         """
diff --git a/trimesh/geometry.py b/trimesh/geometry.py
index 6a9496cc..0c1ae954 100644
--- a/trimesh/geometry.py
+++ b/trimesh/geometry.py
@@ -4,7 +4,7 @@ from . import util
 from .constants import tol, log
 
 try:
-    from scipy.sparse import coo_matrix
+    from scipy.sparse import coo_matrix, identity
 except ImportError:
     log.warning('scipy.sparse.coo_matrix unavailable')
 
@@ -199,6 +199,81 @@ def triangulate_quads(quads):
     return faces
 
 
+def vertex_face_indices(vertex_count,
+                        faces,
+                        **kwargs):
+    """
+    Find vertex face indices from the faces array of vertices
+
+    Parameters
+    -----------
+    vertex_count : int
+      The number of vertices faces refer to
+    faces : (n, 3) int
+      List of vertex indices
+
+    Returns
+    -----------
+    vertex_faces : (vertex_count, ) int
+      Face indices for every vertex
+      Array padded with -1 in each row for all vertices with fewer
+      face indices than the max number of face indices.
+    """
+
+    def sorted_sparse():
+        # use a sparse matrix of which face contains each vertex to
+        # sort the face indices by vertex index
+        # allow cached sparse matrix to be passed
+        if 'sparse' in kwargs:
+            sparse = kwargs['sparse']
+        else:
+            sparse = index_sparse(vertex_count, faces)
+
+        y = identity(len(faces), dtype=int)
+        sorted_faces = sparse.dot(y).nonzero()
+        return sorted_faces
+
+    def sorted_loop():
+        # loop through every vertex, in tests was ~130x slower than
+        # doing this with a sparse matrix
+        sorted_faces = (np.zeros(faces.size), np.zeros(faces.size))
+        for v in np.arange(vertex_count):
+            v_slice = np.arange(vertex_starts[v], vertex_starts[v] + vertex_counts[v])
+            sorted_faces[0][v_slice] = v
+            sorted_faces[1][v_slice] = np.flip(np.where(faces.flatten() == v)[0] // 3, axis=0)
+        return sorted_faces
+
+    # Create 2d array with row for each vertex and
+    # length of max number of faces for a vertex
+    try:
+        vertex_counts = np.bincount(faces.flatten())
+    except TypeError:
+        # casting failed on 32 bit Windows
+        log.error('casting failed!', exc_info=True)
+        # fall back to np.unique (about 35x slower in testing than bincount)
+        _, vertex_counts = np.unique(faces.flatten(), return_counts=True)
+
+    vertex_starts = np.cumsum(vertex_counts)
+    vertex_starts = np.roll(vertex_starts, 1)
+    vertex_starts[0] = 0
+    vertex_faces = np.arange(np.max(vertex_counts)) + vertex_starts[:, None]
+
+    try:
+        sorted_faces = sorted_sparse()
+    except BaseException:
+        log.warning(
+            'unable to generate sparse matrix! Falling back!',
+            exc_info=True)
+        sorted_faces = sorted_loop()
+
+    # Mask each row with number of faces for that vertex and fill with face indices
+    # The rest of the row is filled with -1 to pad to the max length
+    vertex_faces = -(vertex_faces >= (vertex_starts + vertex_counts)[:, None]).astype(int)
+    vertex_faces[vertex_faces == 0] = sorted_faces[1]
+
+    return vertex_faces
+
+
 def mean_vertex_normals(vertex_count,
                         faces,
                         face_normals,
diff --git a/trimesh/grouping.py b/trimesh/grouping.py
index 9306f095..9ed1a18f 100644
--- a/trimesh/grouping.py
+++ b/trimesh/grouping.py
@@ -265,8 +265,9 @@ def unique_ordered(data):
 
 
 def unique_bincount(values,
-                    minlength,
-                    return_inverse=True):
+                    minlength=0,
+                    return_inverse=False,
+                    return_counts=False):
     """
     For arrays of integers find unique values using bin counting.
     Roughly 10x faster for correct input than np.unique
@@ -279,14 +280,19 @@ def unique_bincount(values,
       Maximum value that will occur in values (values.max())
     return_inverse : bool
       If True, return an inverse such that unique[inverse] == values
+    return_counts : bool
+      If True, also return the number of times each unique item appears in values
 
     Returns
     ------------
     unique : (m,) int
       Unique values in original array
-    inverse : (n,) int
+    inverse : (n,) int, optional
       An array such that unique[inverse] == values
       Only returned if return_inverse is True
+    counts : (m,) int, optional
+      An array holding the counts of each unique item in values
+      Only returned if return_counts is True
     """
     values = np.asanyarray(values)
     if len(values.shape) != 1 or values.dtype.kind != 'i':
@@ -299,7 +305,9 @@ def unique_bincount(values,
         # casting failed on 32 bit windows
         log.error('casting failed!', exc_info=True)
         # fall back to numpy unique
-        return np.unique(values, return_inverse=return_inverse)
+        return np.unique(values,
+                         return_inverse=return_inverse,
+                         return_counts=return_counts)
 
     # which bins are occupied at all
     # counts are integers so this works
@@ -308,13 +316,20 @@ def unique_bincount(values,
     # which values are unique
     # indexes correspond to original values
     unique = np.where(unique_bin)[0]
+    ret = (unique,)
 
     if return_inverse:
         # find the inverse to reconstruct original
         inverse = (np.cumsum(unique_bin) - 1)[values]
-        return unique, inverse
+        ret += (inverse,)
 
-    return unique
+    if return_counts:
+        unique_counts = counts[unique]
+        ret += (unique_counts,)
+
+    if len(ret) == 1:
+        return ret[0]
+    return ret
 
 
 def merge_runs(data, digits=None):
diff --git a/trimesh/ray/ray_triangle.py b/trimesh/ray/ray_triangle.py
index f44fb0ee..55762ee1 100644
--- a/trimesh/ray/ray_triangle.py
+++ b/trimesh/ray/ray_triangle.py
@@ -32,43 +32,34 @@ class RayMeshIntersector(object):
                       multiple_hits=True,
                       **kwargs):
         """
-        Find the intersections between the current mesh and an
-        array of rays.
+        Find the intersections between the current mesh and a list of rays.
 
         Parameters
         ------------
-        ray_origins :  (m, 3) float
-          Ray origin points
-        ray_directions : (m, 3) float
-          Ray direction vectors
-        multiple_hits :  bool
-          Consider multiple hits of each ray or not
-        return_locations : bool
-          Return hit locations or not
+        ray_origins:      (m,3) float, ray origin points
+        ray_directions:   (m,3) float, ray direction vectors
+        multiple_hits:    bool, consider multiple hits of each ray or not
+        return_locations: bool, return hit locations or not
 
         Returns
         -----------
-        index_triangle : (h,) int
-          Index of triangles hit
-        index_ray : (h,) int
-          Index of ray that hit triangle
-        locations : (h, 3) float
-          [optional] Position of intersection in space
+        index_triangle: (h,) int,    index of triangles hit
+        index_ray:      (h,) int,    index of ray that hit triangle
+        locations:      (h,3) float, (optional) position of intersection in space
         """
         (index_tri,
          index_ray,
-         locations) = ray_triangle_id(
-             triangles=self.mesh.triangles,
-             ray_origins=ray_origins,
-             ray_directions=ray_directions,
-             tree=self.mesh.triangles_tree,
-             multiple_hits=multiple_hits,
-             triangles_normal=self.mesh.face_normals)
+         locations) = ray_triangle_id(triangles=self.mesh.triangles,
+                                      ray_origins=ray_origins,
+                                      ray_directions=ray_directions,
+                                      tree=self.mesh.triangles_tree,
+                                      multiple_hits=multiple_hits,
+                                      triangles_normal=self.mesh.face_normals)
         if return_locations:
             if len(index_tri) == 0:
                 return index_tri, index_ray, locations
-            unique = grouping.unique_rows(
-                np.column_stack((locations, index_ray)))[0]
+            unique = grouping.unique_rows(np.column_stack((locations,
+                                                           index_ray)))[0]
             return index_tri[unique], index_ray[unique], locations[unique]
         return index_tri, index_ray
 
@@ -116,15 +107,12 @@ class RayMeshIntersector(object):
 
         Parameters
         ------------
-        ray_origins : (m, 3) float
-          Ray origin points
-        ray_directions : (m, 3) float
-          Ray direction vectors
+        ray_origins:      (m,3) float, ray origin points
+        ray_directions:   (m,3) float, ray direction vectors
 
         Returns
         ---------
-        hit : (m,) bool
-          Whether any ray hit any triangle on the mesh
+        hit: boolean, whether any ray hit any triangle on the mesh
         """
         index_tri, index_ray = self.intersects_id(ray_origins,
                                                   ray_directions)
@@ -138,18 +126,15 @@ class RayMeshIntersector(object):
         """
         Check if a mesh contains a list of points, using ray tests.
 
-        If the point is on the surface of the mesh the behavior
-        is undefined.
+        If the point is on the surface of the mesh, behavior is undefined.
 
         Parameters
         ------------
-        points : (n, 3) float
-          Points in space
+        points: (n,3) points in space
 
         Returns
         ---------
-        contains : (n,) bool
-          Whether point is inside mesh or not
+        contains: (n) boolean array, whether point is inside mesh or not
         """
 
         return contains_points(self, points)
@@ -228,27 +213,23 @@ def ray_triangle_id(triangles,
 
     if (len(triangle_candidates) == 0 or
             not valid.any()):
-        # we got no hits so return early with empty array
-        return (np.array([], dtype=np.int64),
-                np.array([], dtype=np.int64),
-                np.array([], dtype=np.float64))
+        return [], [], []
 
     # find the barycentric coordinates of each plane intersection on the
     # triangle candidates
     barycentric = triangles_mod.points_to_barycentric(
         triangle_candidates[valid], location)
 
-    # the plane intersection is inside the triangle if all barycentric
-    # coordinates are between 0.0 and 1.0
+    # the plane intersection is inside the triangle if all barycentric coordinates
+    # are between 0.0 and 1.0
     hit = np.logical_and((barycentric > -tol.zero).all(axis=1),
                          (barycentric < (1 + tol.zero)).all(axis=1))
 
-    # the result index of the triangle is a candidate with a valid
-    # plane intersection and a triangle which contains the plane
-    # intersection point
+    # the result index of the triangle is a candidate with a valid plane intersection and
+    # a triangle which contains the plane intersection point
     index_tri = ray_candidates[valid][hit]
-    # the ray index is a subset with a valid plane intersection and
-    # contained by a triangle
+    # the ray index is a subset with a valid plane intersection and contained
+    # by a triangle
     index_ray = ray_id[valid][hit]
     # locations are already valid plane intersections, just mask by hits
     location = location[hit]
diff --git a/trimesh/version.py b/trimesh/version.py
index df8fe443..7a6bceef 100644
--- a/trimesh/version.py
+++ b/trimesh/version.py
@@ -1,1 +1,1 @@
-__version__ = '2.38.30'
+__version__ = '2.38.29'


Test patch:
diff --git a/tests/test_ray.py b/tests/test_ray.py
index 60bdb53d..94403410 100644
--- a/tests/test_ray.py
+++ b/tests/test_ray.py
@@ -7,8 +7,7 @@ except BaseException:
 class RayTests(g.unittest.TestCase):
 
     def test_rays(self):
-        meshes = [g.get_mesh(**k)
-                  for k in g.data['ray_data']['load_kwargs']]
+        meshes = [g.get_mesh(**k) for k in g.data['ray_data']['load_kwargs']]
         rays = g.data['ray_data']['rays']
         names = [m.metadata['file_name'] for m in meshes]
 
@@ -18,14 +17,14 @@ class RayTests(g.unittest.TestCase):
         for m in meshes:
             name = m.metadata['file_name']
             hit_any.append(m.ray.intersects_any(**rays[name]))
-            hit_loc.append(m.ray.intersects_location(**rays[name])[0])
+            hit_loc.append(len(m.ray.intersects_location(**rays[name])[0]))
             hit_id.append(m.ray.intersects_id(**rays[name]))
         hit_any = g.np.array(hit_any, dtype=g.np.int)
 
         for i in g.trimesh.grouping.group(
                 g.np.unique(names, return_inverse=True)[1]):
             broken = hit_any[i].astype(g.np.int).ptp(axis=0).sum()
-            assert broken == 0
+            self.assertTrue(broken == 0)
 
     def test_rps(self):
         for use_embree in [True, False]:
@@ -40,46 +39,17 @@ class RayTests(g.unittest.TestCase):
             # force ray object to allocate tree before timing it
             # tree = sphere.ray.tree
             tic = [g.time.time()]
-            a = sphere.ray.intersects_id(
-                ray_origins, ray_directions)
+            sphere.ray.intersects_id(ray_origins, ray_directions)
             tic.append(g.time.time())
-            b = sphere.ray.intersects_location(
-                ray_origins, ray_directions)
+            sphere.ray.intersects_location(ray_origins, ray_directions)
             tic.append(g.time.time())
 
-            # make sure ray functions always return numpy arrays
-            assert all(len(i.shape) >= 0 for i in a)
-            assert all(len(i.shape) >= 0 for i in b)
-
             rps = dimension[0] / g.np.diff(tic)
 
             g.log.info('Measured %s rays/second with embree %d',
                        str(rps),
                        use_embree)
 
-    def test_empty(self):
-        """
-        Test queries with no hits
-        """
-        for use_embree in [True, False]:
-            dimension = (100, 3)
-            sphere = g.get_mesh('unit_sphere.STL',
-                                use_embree=use_embree)
-            # should never hit the sphere
-            ray_origins = g.np.random.random(dimension)
-            ray_directions = g.np.tile([0, 1, 0], (dimension[0], 1))
-            ray_origins[:, 2] = -5
-
-            # make sure ray functions always return numpy arrays
-            # these functions return multiple results all of which
-            # should always be a numpy array
-            assert all(len(i.shape) >= 0 for i in
-                       sphere.ray.intersects_id(
-                           ray_origins, ray_directions))
-            assert all(len(i.shape) >= 0 for i in
-                       sphere.ray.intersects_location(
-                           ray_origins, ray_directions))
-
     def test_contains(self):
         scale = 1.5
         for use_embree in [True, False]:
diff --git a/tests/test_util.py b/tests/test_util.py
index 720421e3..d4c267db 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -284,16 +284,22 @@ class UniqueTests(unittest.TestCase):
                 minlength = 10
 
             # try our unique bincount function
-            unique, inverse = g.trimesh.grouping.unique_bincount(
+            unique, inverse, counts = g.trimesh.grouping.unique_bincount(
                 values,
                 minlength=minlength,
-                return_inverse=True)
+                return_inverse=True,
+                return_counts=True)
             # make sure inverse is correct
             assert (unique[inverse] == values).all()
 
+            # make sure that the number of counts matches
+            # the number of unique values
+            assert (len(unique) == len(counts))
+
             # get the truth
-            truth_unique, truth_inverse = np.unique(values,
-                                                    return_inverse=True)
+            truth_unique, truth_inverse, truth_counts = np.unique(values,
+                                                                  return_inverse=True,
+                                                                  return_counts=True)
             # make sure truth is doing what we think
             assert (truth_unique[truth_inverse] == values).all()
 
@@ -303,6 +309,9 @@ class UniqueTests(unittest.TestCase):
             # make sure all values are identical
             assert set(truth_unique) == set(unique)
 
+            # make sure that the truth counts are identical to our counts
+            assert np.all(truth_counts == counts)
+
 
 if __name__ == '__main__':
     trimesh.util.attach_to_log()
diff --git a/tests/test_vertices.py b/tests/test_vertices.py
new file mode 100644
index 00000000..0c4c3da3
--- /dev/null
+++ b/tests/test_vertices.py
@@ -0,0 +1,42 @@
+try:
+    from . import generic as g
+except BaseException:
+    import generic as g
+
+
+class VerticesTest(g.unittest.TestCase):
+
+    def test_vertex_faces(self):
+
+        # One watertight, one not; also various sizes
+        meshes = [g.get_mesh('featuretype.STL'),
+                  g.get_mesh('cycloidal.ply')]
+
+        for m in meshes:
+            # Choose random face indices and make sure that their vertices
+            # contain those face indices
+            rand_faces = g.np.random.randint(low=0, high=len(m.faces), size=100)
+            for f in rand_faces:
+                v_inds = m.faces[f]
+                assert (f in m.vertex_faces[v_inds[0]])
+                assert (f in m.vertex_faces[v_inds[1]])
+                assert (f in m.vertex_faces[v_inds[2]])
+
+            # Choose some random vertices and make sure their face indices
+            # are correct
+            rand_vertices = g.np.random.randint(low=0, high=len(m.vertices), size=100)
+            for v in rand_vertices:
+                v_faces = g.np.flip(g.np.where(m.faces == v)[0], axis=0)
+                assert (g.np.all(v_faces == m.vertex_faces[v][m.vertex_faces[v] >= 0]))
+
+            # Intentionally cause fallback to looping over vertices and make sure we
+            # get the same result
+            loop_vertex_faces = g.trimesh.geometry.vertex_face_indices(len(m.vertices), 
+                                                                       m.faces, 
+                                                                       sparse=None)
+            assert (g.np.all(loop_vertex_faces == m.vertex_faces))
+
+
+if __name__ == '__main__':
+    g.trimesh.util.attach_to_log()
+    g.unittest.main()