generate planar surfaces and slice when slices not provided

reebgraph/reebgraph.py

@@ -5,15 +5,6 @@ import pymeshfix as pf
 import networkx as nx
 
 
-def to_trimesh(mesh):
-    """
-    Converts a mesh into a trimesh format
-    """
-    mesh_fix = pf.MeshFix(mesh)
-    trimesh_mesh = tm.base.Trimesh(vertices=mesh_fix.v, faces=mesh_fix.f)
-    return trimesh_mesh
-
-
 def surf_dists(h0,h1, max_checks=100):
     h0n = h0.compute_normals(point_normals=True, cell_normals=False, auto_orient_normals=True)
 
@@ -37,35 +28,63 @@ def surf_dists(h0,h1, max_checks=100):
     mask = np.isclose(h0n["distances"],0)
     h0n["distances"][mask] = np.nan
     val = np.nanmin(h0n["distances"])
-    print("surf_dist: {}".format(val))
     return val
 
 
+def get_model_xy_grid(mesh, resolution=100):
+    """
+    Get the xy grid used to generate the surfaces for slicing
+    """
+    bounds = mesh.bounds
+    min_x, max_x = bounds[0], bounds[1]
+    min_y, max_y = bounds[2], bounds[3]
+    x = np.linspace(min_x-1e-3, max_x+1e-3, resolution)
+    y = np.linspace(min_y-1e-3, max_y+1e-3, resolution)
+    x, y = np.meshgrid(x, y)
+    return x, y
+
+
+def gen_planar_surf(mesh, z, resolution=100):
+    """
+    Create a slicing surface
+    """
+    x, y = get_model_xy_grid(mesh, resolution)
+    z = np.full([len(x), len(y)], z)
+    surf_points = np.c_[x.reshape(-1), y.reshape(-1), z.reshape(-1)]
+
+    point_cloud = pv.PolyData(surf_points)
+    surf = point_cloud.delaunay_2d()
+    return surf.clean(inplace=False)
+
+
 class ReebGraph:
 
-    def __init__(self, mesh, z_res=10, bydist=False, slices=None):
+    def __init__(self, mesh, planar_zs=[], slices=None):
         """
-        if bydist, z_res is distance between layers
-        otherwise, z_res is number of layers
+        either provide z-vals for planar slicing or provide slices directly
         """
         self.mesh = mesh
-        self.bydist = bydist
-        self.z_res = z_res
+        self.zs = planar_zs
         self.slices=slices
         self.G, self.pv = self._get_reeb()
 
     def _get_reeb(self, min_dist=10):
 
-        tm_mesh = to_trimesh(self.mesh)
+        if self.slices is None and len(self.zs)==0:
+            Exception("Must give non-planar slices or z-heights to planar slice at")
 
+        # if planar surfaces need to be sliced
         if self.slices is None:
-            # TODO do planar slicing
-            if self.bydist:
-                zs = np.arange(tm_mesh.bounds[0,2], tm_mesh.bounds[1,2]+self.z_res, self.z_res)
-            else:
-                zs = np.linspace(tm_mesh.bounds[0,2]+0.01, tm_mesh.bounds[1,2]-0.01, self.z_res)
-            print(zs)
-            self.slices=[None]*len(zs)
+            self.slices = []
+            surf_0 = gen_planar_surf(self.mesh, 0)
+            for i,z in enumerate(self.zs):
+                surf = surf_0.translate([0,0,z], inplace=False)
+                try:
+                    sl = surf.clip_surface(self.mesh, invert=True)
+                    self.slices.append(sl.clean(inplace=False))
+                    print("Slice {} added at z {} of {}".format(i,z,self.zs[-1]))
+                except:
+                    print("No Intersection at height {}".format(z))
 
         G = nx.MultiDiGraph()
         j = 0
@@ -80,7 +99,7 @@ class ReebGraph:
 
                 # sanity check slice size
                 if sl.area <1:
-                    continue # clipping error pyvista
+                    continue # clipping errors in pyvista
 
                 # add node to graph
                 G.add_node(j,
@@ -124,28 +143,28 @@ class ReebGraph:
         pl.show()
 
 
-# THIS IS ONLY HERE FOR TESTING NONPLANAR
-def gen_single_surf(x, y, z):
-    """
-    Generates mesh of a single surface, given xyz data
-    """
-    surf_points = np.c_[x.reshape(-1), y.reshape(-1), z.reshape(-1)]
-    point_cloud = pv.PolyData(surf_points)
-    surf = point_cloud.delaunay_2d()
-    surf = surf.clean(inplace=False)
-    return surf
-def plot_meshes(meshes):
-    """
-    Visualise model meshes, return handle to plotter
-    """
-    plotter = pv.Plotter()
-    for mesh in meshes:
-        plotter.add_mesh(mesh, color=np.random.randint(0, 2, size=3).tolist(), show_edges=True, lighting=False,
-                         opacity=1, line_width=0.1, edge_color=np.random.randint(0, 2, size=3).tolist())
-    plotter.show()
-    return plotter
-
 if __name__ == '__main__':
+    # functions for testing only
+    def gen_single_surf(x, y, z):
+        """
+        Generates mesh of a single surface, given xyz data
+        """
+        surf_points = np.c_[x.reshape(-1), y.reshape(-1), z.reshape(-1)]
+        point_cloud = pv.PolyData(surf_points)
+        surf = point_cloud.delaunay_2d()
+        surf = surf.clean(inplace=False)
+        return surf
+    def plot_meshes(meshes):
+        """
+        Visualise model meshes, return handle to plotter
+        """
+        plotter = pv.Plotter()
+        for mesh in meshes:
+            plotter.add_mesh(mesh, color=np.random.randint(0, 2, size=3).tolist(), show_edges=True, lighting=False,
+                            opacity=1, line_width=0.1, edge_color=np.random.randint(0, 2, size=3).tolist())
+        plotter.show()
+        return plotter
+
     # load mesh
     mesh = pv.read('../horse.stl')
     temp_mesh = pf.MeshFix(mesh)
@@ -156,27 +175,21 @@ if __name__ == '__main__':
 
     # generate non-planar surfs
     bounds = mesh.bounds
-    resolution = 100
-    min_x, max_x = bounds[0], bounds[1]
-    min_y, max_y = bounds[2], bounds[3]
-    min_z, max_z = bounds[4], bounds[5]
-    x = np.linspace(min_x-1e-3, max_x+1e-3, resolution)
-    y = np.linspace(min_y-1e-3, max_y+1e-3, resolution)
-    x, y = np.meshgrid(x, y)
+    x, y = get_model_xy_grid(mesh)
     midx = x[0, x.shape[1] // 2]
     midy = y[y.shape[0] // 2, 0]
     r = ((x[0,0]-x[-1,-1])**2 + (y[0,0]-y[-1,-1])**2)**0.5/1.5
     z_temp = (r**2-(x-midx)**2-(y-midy)**2)**0.5
     z_temp = z_temp - z_temp[y.shape[0]//2,x.shape[1]//2]
-    # z = tu.non_planar_z_arr(x, y, sl.heights, sl.coeffs) if sl.nonplanar else tu.planar_z_arr(x, y, sl.heights)
     surfs = []
-    heights = np.linspace(min_z, max_z+20, 25)
+    heights = np.linspace(mesh.bounds[4], mesh.bounds[5]+20, 25)
     def_surf = gen_single_surf(x, y, z_temp)
     for layer_i in range(len(heights)):
         new_surf = def_surf.copy()
         new_surf.translate([0,0,heights[layer_i]])
         surfs.append(new_surf)
 
+    # generate slices from non-planar surfs
     slices = []
     for layer_i, surf in enumerate(surfs):
         try:
@@ -184,7 +197,8 @@ if __name__ == '__main__':
             slices.append(sl.clean(inplace=False))
             print("Slice", layer_i, "added. At Z", heights[layer_i], "of", heights[-1])
         except:
-            print("Empty Surf")
+            print("No Intersection")
 
-    rg = ReebGraph(mesh, slices=slices)
-    #rg.show()
\ No newline at end of file
+    # gen reeb graphs
+    rg_planar = ReebGraph(mesh, planar_zs=heights)
+    rg_nonplanar = ReebGraph(mesh, slices=slices)
\ No newline at end of file