code: def voxelize_subdivide(mesh,
                         pitch,
                         max_iter=10,
                         edge_factor=2.0):
      
      max_edge = pitch / edge_factor
  
      if max_iter is None:
          longest_edge = np.linalg.norm(mesh.vertices[mesh.edges[:, 0]] -
                                        mesh.vertices[mesh.edges[:, 1]],
                                        axis=1).max()
          max_iter = max(int(np.ceil(np.log2(longest_edge / max_edge))), 0)
  
      # get the same mesh sudivided so every edge is shorter
      # than a factor of our pitch
      v, f = remesh.subdivide_to_size(mesh.vertices,
                                      mesh.faces,
                                      max_edge=max_edge,
                                      max_iter=max_iter)
  
      # convert the vertices to their voxel grid position
      hit = v / pitch
  
      # Provided edge_factor > 1 and max_iter is large enough, this is
      # sufficient to preserve 6-connectivity at the level of voxels.
      hit = np.round(hit).astype(int)
  
      # remove duplicates
      unique, inverse = grouping.unique_rows(hit)
  
      # get the voxel centers in model space
      occupied_index = hit[unique]
  
      origin_index = occupied_index.min(axis=0)
      origin_position = origin_index * pitch
  
      voxels_sparse = (occupied_index - origin_index)
  
      return voxels_sparse, origin_position
docs: Voxelize a surface by subdividing a mesh until every edge is
      shorter than: (pitch / edge_factor)
  
      Parameters
      -----------
      mesh:        Trimesh object
      pitch:       float, side length of a single voxel cube
      max_iter:    int, cap maximum subdivisions or None for no limit.
      edge_factor: float,
  
      Returns
      -----------
      voxels_sparse:   (n,3) int, (m,n,p) indexes of filled cells
      origin_position: (3,) float, position of the voxel
                                   grid origin in space