feat: chunk_mesh (#44)

* wip: mesh chunking (for creating multiple LOD)

* fix: running, though the algorithm is not perfect yet

* fix: off by one error

* fix: use C order

* fix: rounding errors

* fix: allow C or F order

* feat: return grid associations

* fix+perf: mesh.id, switch to all vectors

* test: add basic tests

* fix: getting closer to correctly fixing x direction faces

* fix: able to fix single outliers 6-connected

* refactor: rename to 6 connected

* refactor: add_point and add_triangle to meshobject

* refactor: add const to verts

* fix: support 18 connected

* wip: making progress on dealing with special cases

* fix: a lot closer to working (cube looks ok)

* fix: corners in 18-connected

* fix: precision 3

* fix: sanity check that xaxis is not same as yaxis

* fix: operation should be ceil not rounding

* fix: remove epsilon subtraction

* feat: add concatenate and consolidate to Mesh

* fix: discard not-implemented 3d all different

* refactor: use a single intersection function

* wip: adding support for 26-connected cases

* fix: rounding->ceil on grid size

* fix: consolidate calling wrong name of is_empty

* fix: handle last_face when there are zero verts

* fix: remove debugging statements

* feat: better error messages

* fix: bug in 18-connected code

* perf: faster triangle calculation

* test: add basic chunking test

* fix: correct face winding for 6-connected

* perf: faster face id calculation for 6 connected

* fix: face winding for one case in all_different

* fix: improve windings for 18 connected

Not perfect, but better?

* chore: update test.py into something more useful

* fix: improve windings for close

* feat: automatic conversion to numpy

* fix: safety checks before unsafe array access

* fix: inverted windings

* fix: more inverted windings

* fix: generalize code to handle arbitrary meshes

* fix: two situations for dividing a triangle

* install: fix license metadata and supported python versions

* chore: add namespace to zmesh files

* fix: disallow zero or negative chunk sizes

* fix: handle degenerate intersections

* fix: handle degenerate cases

* wip: attempting to fix the winding logic getting close...

* fix: rejigger below

* fix: incorrect initial offset for faces

* wip: experimenting with a simplified implementation

* fix: fix normals for simple triangle

* fix: check if normals are ccw or cw

* fix: possibly fixed normals

* perf: only compute normal when needed

* refactor: add hat and rename norm to len to Vec3

* refactor: remove premature optimizations

* chore: remove dead code

* fix: remove degenerate faces during consolidate

* refactor: move remove_degenerate_faces into its own function

* test: check for unreferenced vertices

* feat: add remove_unreferenced_vertices

* install: add fastremap to requirements

* docs: add docstring to remove_degenerate_faces

* feat: add removing disconnected vertices to consolidate

* perf: remove unneeded copies

* feat: add merge_close_vertices function

* perf(remove_unreferenced_vertices): much more efficient remapping

* perf(consolidate): remove unnecessary vectorization

* perf: speed up normal remapping

Author: william-silversmith

automated_test.py

@@ -301,4 +301,82 @@ 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_delete_unreference_vertices():
+  vertices = [
+    [0,0,0],
+    [0,1,0],
+    [1,0,0],
+    [5,5,5],
+    [7,7,7],
+    [8,7,7],
+    [7,8,8],
+  ]
+  faces = [[0,1,2],[4,5,6]]
+
+  mesh = zmesh.Mesh(vertices=vertices, faces=faces, normals=None)
+  res = mesh.remove_unreferenced_vertices()
+
+  ans_verts = vertices[:3] + vertices[4:]
+  ans_faces = [[0,1,2],[3,4,5]]
+
+  assert res.vertices.shape == (6,3)
+  assert res.faces.shape == (2,3)
+  assert np.all(res.vertices == np.array(ans_verts))
+  assert np.all(res.faces == np.array(ans_faces))
+
+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
+
+  
+
+
+
+
+
+
+
+
+
 

dev_requirements.txt

@@ -1,4 +1,5 @@
 cython
+fastremap
 pytest
 Jinja2
 tox

setup.cfg

@@ -7,7 +7,7 @@ description_file = README.md
 author = William Silversmith (maintainer), Aleks Zlateski (original author)
 author_email = [email protected]
 home_page = https://github.com/seung-lab/zmesh/
-license = GPL-3.0+
+license = LGPL-3.0-or-later
 classifier =
     Intended Audience :: Developers
     Development Status :: 5 - Production/Stable

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())

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,25 @@ cimport numpy as cnp
 import numpy as np
 from zmesh.mesh import Mesh
 
-cdef extern from "cMesher.hpp":
+cdef extern from "utility.hpp" namespace "zmesh::utility":
   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 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" namespace "zmesh":
   cdef cppclass CMesher[P,L,S]:
     CMesher(vector[float] voxel_res) except +
     void mesh(
@@ -51,6 +66,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 +237,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 +270,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 +372,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,9 @@
 #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
+
+namespace zmesh {
 
 template <typename PositionType, typename LabelType, typename SimplifierType>
 class CMesher {
@@ -49,7 +47,7 @@ class CMesher {
     return marchingcubes_.pack_coords(x,y,z);
   }
 
-  MeshObject get_mesh(
+  zmesh::utility::MeshObject get_mesh(
     LabelType segid, 
     bool generate_normals,
     int simplification_factor,
@@ -60,7 +58,7 @@ class CMesher {
 
     // MC produces no triangles if either
     // none or all voxels were labeled.
-    MeshObject empty_obj;
+    zmesh::utility::MeshObject empty_obj;
     if (marchingcubes_.count(segid) == 0) { 
       return empty_obj;
     }
@@ -77,7 +75,7 @@ class CMesher {
     );
   }
 
-  MeshObject simplify(      
+  zmesh::utility::MeshObject simplify(      
       const std::vector< zi::vl::vec< PositionType, 3> >& triangles,
       bool generate_normals,
       int simplification_factor,
@@ -89,7 +87,7 @@ class CMesher {
       bool transpose = true 
     ) {
 
-    MeshObject obj;
+    zmesh::utility::MeshObject obj;
 
     zi::mesh::int_mesh<PositionType, LabelType> im;
     im.add(triangles);
@@ -177,7 +175,7 @@ class CMesher {
     return obj;
   }
 
-  MeshObject simplify_points(
+  zmesh::utility::MeshObject simplify_points(
     const uint64_t* points,
     const size_t Nv,
     bool generate_normals,
@@ -219,5 +217,6 @@ class CMesher {
   }
 };
 
+};
 
 #endif