Repository: mikedh/trimesh
Pull request #260 (2018-11-24 21:12:50)

Fix indices_to_points to make it consistent to VoxelMesh.as_boxes
I'm not so sure, but usually voxel point coordinate should be calculated as:

  (index - 0.5) * pitch + origin = index * pitch + origin - pitch / 2.0

In VoxelMesh.as_boxes this formulation is used, but the matrix_to_points function does not.


Patch:
diff --git a/trimesh/nsphere.py b/trimesh/nsphere.py
index 06fc09af..c469b6fc 100644
--- a/trimesh/nsphere.py
+++ b/trimesh/nsphere.py
@@ -17,9 +17,6 @@ try:
 except ImportError:
     log.warning('No scipy!')
 
-# if we get a MemoryError change this to True
-# so we don't keep trying tile operations
-_LOWMEM = False
 
 def minimum_nsphere(obj):
     """
@@ -70,35 +67,25 @@ def minimum_nsphere(obj):
     # , qhull_options='QbB Pp')
     voronoi = spatial.Voronoi(points, furthest_site=True)
 
-    # if we've gotten a MemoryError before don't try tiling
-    global _LOWMEM
-
     # find the maximum radius^2 point for each of the voronoi vertices
     # this is worst case quite expensive, but we have used quick convex
     # hull methods to reduce n for this operation
     # we are doing comparisons on the radius squared then rooting once
-    if _LOWMEM or (len(points) * len(voronoi.vertices)) > 1e7:
+
+    if len(points) * len(voronoi.vertices) > 1e7:
         # if we have a bajillion points loop
         radii_2 = np.array([((points - v)**2).sum(axis=1).max()
                             for v in voronoi.vertices])
     else:
         # otherwise tiling is massively faster
-        try:
-            dim = points.shape[1]
-            v_tile = np.tile(voronoi.vertices, len(points)).reshape((-1, dim))
-            # tile points per voronoi vertex
-            p_tile = np.tile(points, (len(voronoi.vertices), 1)).reshape((-1, dim))
-            # find the maximum radius of points for each voronoi vertex
-            radii_2 = ((p_tile - v_tile)**2).sum(axis=1).reshape(
-                (len(voronoi.vertices), -1)).max(axis=1)
-        except MemoryError:
-            # don't try tiling again
-            _LOWMEM = True
-            # fall back to the list comprehension
-            radii_2 = np.array([((points - v)**2).sum(axis=1).max()
-                                for v in voronoi.vertices])
-            # log the MemoryError
-            log.warning('MemoryError: falling back to slower check!')
+        # what dimension are the points
+        dim = points.shape[1]
+        v_tile = np.tile(voronoi.vertices, len(points)).reshape((-1, dim))
+        # tile points per voronoi vertex
+        p_tile = np.tile(points, (len(voronoi.vertices), 1)).reshape((-1, dim))
+        # find the maximum radius of points for each voronoi vertex
+        radii_2 = ((p_tile - v_tile)**2).sum(axis=1).reshape(
+            (len(voronoi.vertices), -1)).max(axis=1)
 
     # we want the smallest sphere, so we take the min of the radii options
     radii_idx = radii_2.argmin()
diff --git a/trimesh/version.py b/trimesh/version.py
index ef754220..d2190d98 100644
--- a/trimesh/version.py
+++ b/trimesh/version.py
@@ -1,1 +1,1 @@
-__version__ = '2.35.35'
+__version__ = '2.35.34'
diff --git a/trimesh/voxel.py b/trimesh/voxel.py
index 966c48b3..c03121f8 100644
--- a/trimesh/voxel.py
+++ b/trimesh/voxel.py
@@ -14,7 +14,7 @@ from . import grouping
 from .constants import log, log_time
 
 
-class VoxelBase(object):
+class Voxel(object):
 
     def __init__(self, *args, **kwargs):
         self._data = caching.DataStore()
@@ -109,12 +109,10 @@ class VoxelBase(object):
         ---------
         index: (3,) int tuple, index in self.matrix
         """
-        point = np.asanyarray(point)
-        if point.shape != (3,):
-            raise ValueError('to_index requires a single point')
-        index = np.round((point - self.origin) /
-                         self.pitch).astype(int)
-        index = tuple(index)
+        indices = points_to_indices(points=[point],
+                                    pitch=self.pitch,
+                                    origin=self.origin)
+        index = tuple(indices[0])
         return index
 
     def is_filled(self, point):
@@ -138,49 +136,7 @@ class VoxelBase(object):
         return is_filled
 
 
-class Voxel(VoxelBase):
-
-    def __init__(self, matrix, pitch, origin):
-        super(Voxel, self).__init__()
-
-        self._data['matrix'] = matrix
-        self._data['pitch'] = pitch
-        self._data['origin'] = origin
-
-    @property
-    def origin(self):
-        return self._data['origin']
-
-    @property
-    def matrix(self):
-        return self._data['matrix']
-
-    def as_boxes(self):
-        """
-        A rough Trimesh representation of the voxels with a box
-        for each filled voxel.
-
-        Returns
-        ---------
-        mesh: Trimesh object made up of one box per filled cell.
-        """
-        centers = matrix_to_points(
-            matrix=self._data['matrix'],
-            pitch=self._data['pitch'],
-            origin=self._data['origin'],
-        )
-        mesh = multibox(centers=centers, pitch=self.pitch)
-        return mesh
-
-    def show(self, *args, **kwargs):
-        """
-        Convert the current set of voxels into a trimesh for visualization
-        and show that via its built- in preview method.
-        """
-        return self.as_boxes().show(*args, **kwargs)
-
-
-class VoxelMesh(VoxelBase):
+class VoxelMesh(Voxel):
 
     def __init__(self,
                  mesh,
@@ -316,8 +272,9 @@ class VoxelMesh(VoxelBase):
         else:
             filled = self.sparse_surface
         # center points of voxels
-        centers = (filled * self.pitch).astype(np.float64)
-        centers += self.origin - (self.pitch / 2.0)
+        centers = indices_to_points(indices=filled,
+                                    pitch=self.pitch,
+                                    origin=self.origin)
         mesh = multibox(centers=centers, pitch=self.pitch)
         return mesh
 
@@ -603,6 +560,62 @@ def fill_voxelization(occupied):
     return filled
 
 
+def points_to_indices(points, pitch, origin):
+    """
+    Convert center points of an (n,m,p) matrix into its indices.
+
+    Parameters
+    ----------
+    points: (q, 3) float, center points of voxel matrix (n,m,p)
+    pitch: float, what pitch was the voxel matrix computed with
+    origin: (3,) float, what is the origin of the voxel matrix
+
+    Returns
+    ----------
+    indices: (q, 3) int, list of indices
+    """
+    points = np.asanyarray(points, dtype=np.float64)
+    origin = np.asanyarray(origin, dtype=np.float64)
+    pitch = float(pitch)
+
+    if points.shape != (points.shape[0], 3):
+        raise ValueError('shape of points must be (q, 3)')
+
+    if origin.shape != (3,):
+        raise ValueError('shape of origin must be (3,)')
+
+    indices = np.round((points - origin) / pitch + 0.5).astype(int)
+    return indices
+
+
+def indices_to_points(indices, pitch, origin):
+    """
+    Convert indices of an (n,m,p) matrix into a set of voxel center points.
+
+    Parameters
+    ----------
+    indices: (q, 3) int, index of voxel matrix (n,m,p)
+    pitch: float, what pitch was the voxel matrix computed with
+    origin: (3,) float, what is the origin of the voxel matrix
+
+    Returns
+    ----------
+    points: (q, 3) float, list of points
+    """
+    indices = np.asanyarray(indices, dtype=np.float64)
+    origin = np.asanyarray(origin, dtype=np.float64)
+    pitch = float(pitch)
+
+    if indices.shape != (indices.shape[0], 3):
+        raise ValueError('shape of indices must be (q, 3)')
+
+    if origin.shape != (3,):
+        raise ValueError('shape of origin must be (3,)')
+
+    points = (indices - 0.5) * pitch + origin
+    return points
+
+
 def matrix_to_points(matrix, pitch, origin):
     """
     Convert an (n,m,p) matrix into a set of points for each voxel center.
@@ -617,7 +630,8 @@ def matrix_to_points(matrix, pitch, origin):
     ----------
     points: (q, 3) list of points
     """
-    points = np.column_stack(np.nonzero(matrix)) * pitch + origin
+    indices = np.column_stack(np.nonzero(matrix))
+    points = indices_to_points(indices=indices, pitch=pitch, origin=origin)
     return points
 
 


Test patch:
diff --git a/tests/test_voxel.py b/tests/test_voxel.py
index e27083d4..5fe54e69 100644
--- a/tests/test_voxel.py
+++ b/tests/test_voxel.py
@@ -101,6 +101,33 @@ class VoxelTest(g.unittest.TestCase):
                                                 radius=2,
                                                 fill=True)
 
+    def test_points_to_from_indices(self):
+        points = [[0.05, 0.55, 0.35]]
+        origin = [0, 0, 0]
+        pitch = 0.1
+        indices = [[1, 6, 4]]
+
+        # points -> indices
+        indices2 = g.trimesh.voxel.points_to_indices(
+            points=points, origin=origin, pitch=pitch
+        )
+        g.np.testing.assert_allclose(indices, indices2, atol=0, rtol=0)
+
+        # indices -> points
+        points2 = g.trimesh.voxel.indices_to_points(
+            indices=indices, origin=origin, pitch=pitch
+        )
+        g.np.testing.assert_allclose(points, points2)
+
+        # indices -> points -> indices (this must be consistent)
+        points2 = g.trimesh.voxel.indices_to_points(
+            indices=indices, origin=origin, pitch=pitch
+        )
+        indices2 = g.trimesh.voxel.points_to_indices(
+            points=points2, origin=origin, pitch=pitch
+        )
+        g.np.testing.assert_allclose(indices, indices2, atol=0, rtol=0)
+
 
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()