feat: chunk_mesh

automated_test.py

@@ -296,4 +296,59 @@ def test_min_error_skip(reduction_factor):
     assert abs(factor - reduction_factor) < 1
 
 
+def test_chunk_shape():
+  labels = np.zeros( (10, 10, 10), dtype=np.uint32)
+  labels[1:-1, 1:-1, 1:-1] = 1
+
+  mesher = zmesh.Mesher( (1, 1, 1) )
+  mesher.mesh(labels)
+  mesh = mesher.get_mesh(1, normals=False, simplification_factor=0, max_simplification_error=100)
+
+  meshes = zmesh.chunk_mesh(
+    mesh,
+    [5.,5.,5.],
+  )
+
+  assert len(meshes) == 8
+  assert not any([
+    m.is_empty() for m in meshes.values()
+  ])
+
+def test_chunk_mesh_triangle():
+  vertices = [
+    [0,0,0],
+    [0,1,0],
+    [1,0,0],
+  ]
+  faces = [[0,1,2]]
+
+  mesh = zmesh.Mesh(vertices=vertices, faces=faces, normals=None)
+
+  meshes = zmesh.chunk_mesh(mesh, [.5,.5,.5])
+
+  meshes = [ m for m in meshes.values() if not m.is_empty() ]
+
+  assert len(meshes) == 3
+
+  m = zmesh.Mesh.concatenate(*meshes).consolidate()
+
+  assert m.vertices.shape[0] == 6
+
+  assert [0,0,0] in m.vertices
+  assert [1,0,0] in m.vertices
+  assert [0,1,0] in m.vertices
+  assert [0,0.5,0] in m.vertices
+  assert [0.5,0,0] in m.vertices
+  assert [0.5,0.5,0] in m.vertices
+
+
+
+
+
+
+
+
+
+
+
 

test.py

@@ -1,18 +1,19 @@
 import numpy as np
 import zmesh
 
-labels = np.zeros( (50, 50, 50), dtype=np.uint32)
+labels = np.zeros( (20, 20, 20), dtype=np.uint32)
 labels[1:-1, 1:-1, 1:-1] = 1
 
-mesher = zmesh.Mesher( (3.141, 1, 1) )
+mesher = zmesh.Mesher( (1, 1, 1) )
 mesher.mesh(labels)
-mesh = mesher.get_mesh(1, normals=False, simplification_factor=100, max_simplification_error=100)
+mesh = mesher.get_mesh(1, normals=False, simplification_factor=0, max_simplification_error=100)
 
-print(mesh)
 
-with open('wow.obj', 'bw') as f:
-  f.write(mesh.to_obj())
-
-with open('wow.ply', 'bw') as f:
-  f.write(mesh.to_ply())
+meshes = zmesh.chunk_mesh(
+  mesh,
+  [10.,10.,10.],
+)
+m = zmesh.Mesh.concatenate(*list(meshes.values()))
 
+with open('cube.obj', 'bw') as f:
+  f.write(m.to_obj())

zi_lib/zi/mesh/quadratic_simplifier.hpp

@@ -468,7 +468,7 @@ class simplifier : non_copyable
 
         // double no_faces = static_cast< double >( mesh_.face_count() );
 
-        std::size_t bad = 0;
+        // std::size_t bad = 0;
         while (heap_.size())
         {
             if (((mesh_.face_count() <= target_faces) &&
@@ -477,10 +477,13 @@ class simplifier : non_copyable
             {
                 break;
             }
-            if (!iterate())
-            {
-                ++bad;
-            }
+
+            iterate();
+
+            // if (!iterate())
+            // {
+            //     ++bad;
+            // }
         }
 
         // generate_normals();

zmesh/_zmesh.pyx

@@ -2,6 +2,8 @@
 # distutils: language = c++
 import cython
 
+from typing import List, Tuple, Sequence
+
 from libc.stdint cimport uint64_t, uint32_t, uint16_t, uint8_t
 from libcpp.vector cimport vector
 from libcpp cimport bool
@@ -14,12 +16,19 @@ cimport numpy as cnp
 import numpy as np
 from zmesh.mesh import Mesh
 
-cdef extern from "cMesher.hpp":
+cdef extern from "utility.hpp":
   cdef struct MeshObject:
     vector[float] points
     vector[float] normals
     vector[unsigned int] faces
 
+  cdef vector[MeshObject] chunk_mesh_accelerated(
+    float* vertices, uint64_t num_vertices,
+    unsigned int* faces, uint64_t num_faces,
+    float cx, float cy, float cz
+  ) except +
+
+cdef extern from "cMesher.hpp":
   cdef cppclass CMesher[P,L,S]:
     CMesher(vector[float] voxel_res) except +
     void mesh(
@@ -51,6 +60,78 @@ cdef extern from "cMesher.hpp":
     void clear()
     P pack_coords(P x, P y, P z)
 
+def chunk_mesh(
+  mesh:Mesh,
+  chunk_size:Sequence[float],
+) -> Dict[Tuple[int,int,int], Mesh]:
+
+  vert_order = 'C' if mesh.vertices.flags.c_contiguous else 'F'
+  face_order = 'C' if mesh.faces.flags.c_contiguous else 'F'
+
+  cdef cnp.ndarray[float] vertices = mesh.vertices.reshape([mesh.vertices.size], order=vert_order)
+  cdef cnp.ndarray[unsigned int] faces = mesh.faces.reshape([mesh.faces.size], order=face_order)
+
+  if vertices.size % 3 != 0:
+    raise ValueError(f"Invalid vertex array. Must be a multiple of 3. Got: {vertices.size}")
+
+  if faces.size % 3 != 0:
+    raise ValueError(f"Invalid faces array. Must be a multiple of 3. Got: {faces.size}")
+
+  cdef vector[MeshObject] objs = chunk_mesh_accelerated(
+    <float*>&vertices[0], mesh.vertices.shape[0],
+    <unsigned int*>&faces[0], mesh.faces.shape[0],
+    chunk_size[0], chunk_size[1], chunk_size[2]  
+  )
+
+  def norm(msh):
+    points = np.array(msh['points'], dtype=np.float32)
+    Nv = points.size // 3
+    Nf = len(msh['faces']) // 3
+
+    points = points.reshape(Nv, 3)
+    faces = np.array(msh['faces'], dtype=np.uint32).reshape(Nf, 3)
+    m = Mesh(points, faces, None)
+    if hasattr(msh, 'id'):
+      m.id = msh.id
+    return m
+
+  minpt = np.min(mesh.vertices, axis=0)
+  maxpt = np.max(mesh.vertices, axis=0)
+
+  grid_size = np.ceil((maxpt - minpt) / np.array(chunk_size)).astype(int)
+  grid_size = np.maximum(grid_size, [1,1,1])
+
+  chunked_meshes = {}
+  i = 0
+  for gz in range(grid_size[2]):
+    for gy in range(grid_size[1]):
+      for gx in range(grid_size[0]):
+        chunked_meshes[(gx,gy,gz)] = norm(objs[i])
+        i += 1
+
+  return chunked_meshes
+
+def _normalize_mesh(mesh, voxel_centered, physical, resolution):
+  """Convert a MeshObject into a  zmesh.Mesh."""
+  points = np.array(mesh['points'], dtype=np.float32)
+  Nv = points.size // 3
+  Nf = len(mesh['faces']) // 3
+
+  points = points.reshape(Nv, 3)
+  if not physical:
+    points *= resolution
+
+  if voxel_centered:
+    points += resolution
+  points /= 2.0
+  faces = np.array(mesh['faces'], dtype=np.uint32).reshape(Nf, 3)
+
+  normals = None
+  if mesh['normals']:
+    normals = np.array(mesh['normals'], dtype=np.float32).reshape(Nv, 3)
+
+  return Mesh(points, faces, normals)
+
 class Mesher:
   """
   Represents a meshed volume. 
@@ -150,7 +231,7 @@ class Mesher:
       transpose=True
     )
 
-    return self._normalize_simplified_mesh(mesh, voxel_centered, physical=True)
+    return _normalize_mesh(mesh, voxel_centered, physical=True, resolution=self.voxel_res)
 
   @cython.binding(True)
   def get(
@@ -183,30 +264,10 @@ class Mesher:
       transpose=False
     )
 
-    mesh = self._normalize_simplified_mesh(mesh, voxel_centered, physical=True)
+    mesh = _normalize_mesh(mesh, voxel_centered, physical=True, resolution=self.voxel_res)
     mesh.id = int(label)
     return mesh
 
-  def _normalize_simplified_mesh(self, mesh, voxel_centered, physical):
-    points = np.array(mesh['points'], dtype=np.float32)
-    Nv = points.size // 3
-    Nf = len(mesh['faces']) // 3
-
-    points = points.reshape(Nv, 3)
-    if not physical:
-      points *= self.voxel_res
-
-    if voxel_centered:
-      points += self.voxel_res
-    points /= 2.0
-    faces = np.array(mesh['faces'], dtype=np.uint32).reshape(Nf, 3)
-
-    normals = None
-    if mesh['normals']:
-      normals = np.array(mesh['normals'], dtype=np.float32).reshape(Nv, 3)
-
-    return Mesh(points, faces, normals)
-
   @cython.binding(True)
   def compute_normals(self, mesh):
     """
@@ -305,7 +366,7 @@ class Mesher:
       for i in range(len(result.points)):
         result.points[i] += min_vertex
 
-    return self._normalize_simplified_mesh(result, voxel_centered, physical=False)
+    return _normalize_mesh(result, voxel_centered, physical=False, resolution=self.voxel_res)
 
   def clear(self):
     self._mesher.clear()

zmesh/cMesher.hpp

@@ -7,11 +7,7 @@
 #include <zi/mesh/quadratic_simplifier.hpp>
 #include <zi/vl/vec.hpp>
 
-struct MeshObject {
-  std::vector<float> points;
-  std::vector<float> normals;
-  std::vector<unsigned int> faces;
-};
+#include "utility.hpp" // includes MeshObject def'n
 
 template <typename PositionType, typename LabelType, typename SimplifierType>
 class CMesher {

zmesh/mesh.py

@@ -14,9 +14,14 @@ class Mesh:
 
   """
   def __init__(self, vertices, faces, normals, id=None):
-    self.vertices = vertices
-    self.faces = faces
-    self.normals = normals
+    self.vertices = np.asarray(vertices, dtype=np.float32)
+    self.faces = np.asarray(faces, dtype=np.uint32)
+
+    if normals is None:
+      self.normals = normals
+    else:
+      self.normals = np.asarray(normals, dtype=np.float32)
+
     self.id = id
 
   def __len__(self):
@@ -45,6 +50,9 @@ def __repr__(self):
       (None if self.normals is None else self.normals.shape[0])
     )
 
+  def is_empty(self):
+    return self.faces.size == 0 or self.vertices.size == 0
+
   @property
   def nbytes(self) -> int:
     nbytes = self.vertices.nbytes if self.vertices is not None else 0
@@ -58,15 +66,46 @@ def clone(self):
     else:
       return Mesh(np.copy(self.vertices), np.copy(self.faces), np.copy(self.normals))
 
-  def triangles(self):
-    Nf = self.faces.shape[0]
-    tris = np.zeros( (Nf, 3, 3), dtype=np.float32, order='C' ) # triangle, vertices, (x,y,z)
+  def triangles(self) -> np.ndarray:
+    """Returns vertex triples representing triangluar faces."""
+    return self.vertices[self.faces]
+
+  @classmethod
+  def concatenate(cls, *meshes, id=None):
+    vertex_ct = np.zeros(len(meshes) + 1, np.uint32)
+    vertex_ct[1:] = np.cumsum([ len(mesh) for mesh in meshes ])
+
+    vertices = np.concatenate([ mesh.vertices for mesh in meshes ])
+
+    faces = np.concatenate([ 
+      mesh.faces + vertex_ct[i] for i, mesh in enumerate(meshes) 
+    ])
+
+    # normals = np.concatenate([ mesh.normals for mesh in meshes ])
+
+    return Mesh(vertices, faces, None, id=id)
+
+  def consolidate(self):
+    """Remove duplicate vertices and faces. Returns a new mesh object."""
+    if self.is_empty():
+      return Mesh([], [], normals=None)
+
+    vertices = self.vertices
+    faces = self.faces
+    normals = self.normals
 
-    for i in range(Nf):
-      for j in range(3):
-        tris[i,j,:] = self.vertices[ self.faces[i,j] ]
+    eff_verts, uniq_idx, idx_representative = np.unique(
+      vertices, axis=0, return_index=True, return_inverse=True
+    )
+
+    face_vector_map = np.vectorize(lambda x: idx_representative[x])
+    eff_faces = face_vector_map(faces)
+    eff_faces = np.unique(eff_faces, axis=0)
 
-    return tris
+    # normal_vector_map = np.vectorize(lambda idx: normals[idx])
+    # eff_normals = normal_vector_map(uniq_idx)
+
+    return Mesh(eff_verts, eff_faces, None, id=self.id)
 
   @classmethod
   def from_precomputed(self, binary):
@@ -231,7 +270,7 @@ def viewer(self):
     render_window_interactor.SetRenderWindow(render_window)
 
     render_window.SetSize(1024, 1024)
-
+      
     renderer.AddActor(actor)
     renderer.SetBackground(0.1, 0.2, 0.3)  # Background color