Sparse Voxel Grid (#937)

* Add a utility function which computes a sparse set of epsilon-sized grid cells surrounding an input surface. The function takes a point on the surface as input and finds all those epsilon grid cells that intersect the surface using a breadth first search over the grid space.

* Add a utility function which computes a sparse set of epsilon-sized grid cells surrounding an input surface. The function takes a point on the surface as input and finds all those epsilon grid cells that intersect the surface using a breadth first search over the grid space.


Former-commit-id: 4f030acb671a6d1499348df76a097da4a35aaaab

include/igl/copyleft/marching_cubes.h

@@ -65,7 +65,7 @@ namespace igl
 
     // marching_cubes( values, points, indices, vertices, faces )
     //
-    // Perform marching cubes reconstruction on the grid cells defined by indices.
+    // Perform marching cubes reconstruction on the sparse grid cells defined by (indices, points).
     // The indices parameter is an nx8 dense array of index values into the points and values arrays.
     // Each row of indices represents a cube for which to generate vertices and faces over.
     //

include/igl/sparse_voxel_grid.cpp

@@ -0,0 +1,140 @@
+#include "sparse_voxel_grid.h"
+
+#include <unordered_map>
+#include <array>
+#include <vector>
+
+
+template <typename DerivedS, typename DerivedP0, typename DerivedV, typename DerivedI>
+IGL_INLINE void igl::sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
+                                       const std::function<typename DerivedS::Scalar(const DerivedP0&)>& scalarFunc,
+                                       const double eps,
+                                       Eigen::PlainObjectBase<DerivedS>& CS,
+                                       Eigen::PlainObjectBase<DerivedV>& CV,
+                                       Eigen::PlainObjectBase<DerivedI>& CI,
+                                       int expected_number_of_cubes)
+{
+  typedef typename DerivedV::Scalar ScalarV;
+  typedef typename DerivedS::Scalar ScalarS;
+  typedef typename DerivedI::Scalar ScalarI;
+  typedef Eigen::Matrix<ScalarV, 1, 3> VertexRowVector;
+  typedef Eigen::Matrix<ScalarI, 1, 8> IndexRowVector;
+
+
+  struct IndexRowVectorHash {
+    std::size_t operator()(const Eigen::RowVector3i& key) const {
+      std::size_t seed = 0;
+      std::hash<int> hasher;
+      for (int i = 0; i < 3; i++) {
+        seed ^= hasher(key[i]) + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine
+      }
+      return seed;
+    }
+  };
+
+  auto sgn = [](ScalarS val) -> int {
+    return (ScalarS(0) < val) - (val < ScalarS(0));
+  };
+
+  ScalarV half_eps = 0.5 * eps;
+
+  std::vector<IndexRowVector> CI_vector(expected_number_of_cubes);
+  std::vector<VertexRowVector> CV_vector(8*expected_number_of_cubes);
+  std::vector<ScalarS> CS_vector(8*expected_number_of_cubes);
+
+  // Track visisted neighbors
+  std::unordered_map<Eigen::RowVector3i, int, IndexRowVectorHash> visited(6*expected_number_of_cubes);
+
+  // BFS Queue
+  std::vector<Eigen::RowVector3i> queue(expected_number_of_cubes*8);
+  queue.push_back(Eigen::RowVector3i(0, 0, 0));
+  while (queue.size() > 0)
+  {
+    Eigen::RowVector3i pi = queue.back();
+    queue.pop_back();
+
+    VertexRowVector ctr = p0 + eps*pi.cast<ScalarV>(); // R^3 center of this cube
+
+    // X, Y, Z basis vectors, and array of neighbor offsets used to construct cubes
+    const Eigen::RowVector3i bx(1, 0, 0), by(0, 1, 0), bz(0, 0, -1);
+    const std::array<Eigen::RowVector3i, 6> neighbors = {
+      bx, -bx, by, -by, bz, -bz
+    };
+
+    // Compute the position of the cube corners and the scalar values at those corners
+    std::array<VertexRowVector, 8> cubeCorners = {
+      ctr+half_eps*(bx+by+bz).cast<ScalarV>(), ctr+half_eps*(bx+by-bz).cast<ScalarV>(), ctr+half_eps*(-bx+by-bz).cast<ScalarV>(), ctr+half_eps*(-bx+by+bz).cast<ScalarV>(),
+      ctr+half_eps*(bx-by+bz).cast<ScalarV>(), ctr+half_eps*(bx-by-bz).cast<ScalarV>(), ctr+half_eps*(-bx-by-bz).cast<ScalarV>(), ctr+half_eps*(-bx-by+bz).cast<ScalarV>()
+    };
+    std::array<ScalarS, 8> cubeScalars;
+    for (int i = 0; i < 8; i++) { cubeScalars[i] = scalarFunc(cubeCorners[i]); }
+
+    // If this cube doesn't intersect the surface, disregard it
+    bool validCube = false;
+    int sign = sgn(cubeScalars[0]);
+    for (int i = 1; i < 8; i++) {
+      if (sign != sgn(cubeScalars[i])) {
+        validCube = true;
+        break;
+      }
+    }
+    if (!validCube) {
+      continue;
+    }
+
+    // Add the cube vertices and indices to the output arrays if they are not there already
+    IndexRowVector cube;
+    uint8_t vertexAlreadyAdded = 0; // This is a bimask. If a bit is 1, it has been visited already by the BFS
+    constexpr std::array<uint8_t, 6> zv = {
+      (1 << 0) | (1 << 1) | (1 << 4) | (1 << 5),
+      (1 << 2) | (1 << 3) | (1 << 6) | (1 << 7),
+      (1 << 0) | (1 << 1) | (1 << 2) | (1 << 3),
+      (1 << 4) | (1 << 5) | (1 << 6) | (1 << 7),
+      (1 << 0) | (1 << 3) | (1 << 4) | (1 << 7),
+      (1 << 1) | (1 << 2) | (1 << 5) | (1 << 6), };
+    constexpr std::array<std::array<int, 4>, 6> zvv {{
+        {{0, 1, 4, 5}}, {{3, 2, 7, 6}}, {{0, 1, 2, 3}},
+        {{4, 5, 6, 7}}, {{0, 3, 4, 7}}, {{1, 2, 5, 6}} }};
+
+    for (int n = 0; n < 6; n++) { // For each neighbor, check the hash table to see if its been added before
+      Eigen::RowVector3i nkey = pi + neighbors[n];
+      auto nbr = visited.find(nkey);
+      if (nbr != visited.end()) { // We've already visited this neighbor, use references to its vertices instead of duplicating them
+        vertexAlreadyAdded |= zv[n];
+        for (int i = 0; i < 4; i++) { cube[zvv[n][i]] = CI_vector[nbr->second][zvv[n % 2 == 0 ? n + 1 : n - 1][i]]; }
+      } else {
+        queue.push_back(nkey); // Otherwise, we have not visited the neighbor, put it in the BFS queue
+      }
+    }
+
+    for (int i = 0; i < 8; i++) { // Add new, non-visited,2 vertices to the arrays
+      if (0 == ((1 << i) & vertexAlreadyAdded)) {
+        cube[i] = CS_vector.size();
+        CV_vector.push_back(cubeCorners[i]);
+        CS_vector.push_back(cubeScalars[i]);
+      }
+    }
+
+    visited[pi] = CI_vector.size();
+    CI_vector.push_back(cube);
+  }
+
+  CV.conservativeResize(CV_vector.size(), 3);
+  CS.conservativeResize(CS_vector.size(), 1);
+  CI.conservativeResize(CI_vector.size(), 8);
+  // If you pass in column-major matrices, this is going to be slooooowwwww
+  for (int i = 0; i < CV_vector.size(); i++) {
+    CV.row(i) = CV_vector[i];
+  }
+  for (int i = 0; i < CS_vector.size(); i++) {
+    CS[i] = CS_vector[i];
+  }
+  for (int i = 0; i < CI_vector.size(); i++) {
+    CI.row(i) = CI_vector[i];
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::sparse_voxel_grid<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, std::function<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&)> const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, int);
+#endif

include/igl/sparse_voxel_grid.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Francis Williams <francis@fwilliams.info>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SPARSE_VOXEL_GRID_H
+#define IGL_SPARSE_VOXEL_GRID_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+namespace igl {
+
+  // sparse_voxel_grid( p0, scalarFunc, eps, CV, CS, CI )
+  //
+  // Given a point, p0, on an isosurface, construct a shell of epsilon sized cubes surrounding the surface.
+  // These cubes can be used as the input to marching cubes.
+  //
+  // Input:
+  //  p0  A 3D point on the isosurface surface defined by scalarFunc(x) = 0
+  //  scalarFunc  A scalar function from R^3 to R -- points which map to 0 lie
+  //              on the surface, points which are negative lie inside the surface,
+  //              and points which are positive lie outside the surface
+  //  eps  The edge length of the cubes surrounding the surface
+  //  expected_number_of_cubes  This pre-allocates internal data structures to speed things up
+  // Output:
+  //   CS  #cube-vertices by 1 list of scalar values at the cube vertices
+  //   CV  #cube-vertices by 3 list of cube vertex positions
+  //   CI  #number of cubes by 8 list of indexes into CS and CV. Each row represents a cube
+  //
+  template <typename DerivedS, typename DerivedP0, typename DerivedV, typename DerivedI>
+  IGL_INLINE void sparse_voxel_grid(const Eigen::MatrixBase<DerivedP0>& p0,
+                                    const std::function<typename DerivedS::Scalar(const DerivedP0&)>& scalarFunc,
+                                    const double eps,
+                                    Eigen::PlainObjectBase<DerivedS>& CS,
+                                    Eigen::PlainObjectBase<DerivedV>& CV,
+                                    Eigen::PlainObjectBase<DerivedI>& CI,
+                                    int expected_number_of_cubes=1024);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#    include "sparse_voxel_grid.cpp"
+#endif
+
+#endif // IGL_SPARSE_VOXEL_GRID_H

tutorial/715_MeshImplicitFunction/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)

tutorial/715_MeshImplicitFunction/main.cpp

@@ -0,0 +1,50 @@
+#include <igl/copyleft/marching_cubes.h>
+#include <igl/sparse_voxel_grid.h>
+#include <igl/opengl/glfw/Viewer.h>
+
+#include <Eigen/Core>
+#include <iostream>
+
+#include "tutorial_shared_path.h"
+
+int main(int argc, char * argv[])
+{
+  // An implicit function which is zero on the surface of a sphere centered at the origin with radius 1
+  // This function is negative inside the surface and positive outside the surface
+  std::function<double(const Eigen::RowVector3d&)> scalar_func = [](const Eigen::RowVector3d& pt) -> double {
+    return pt.norm() - 1.0;
+  };
+
+  // We know that the point (0, 0, 1) lies on the implicit surface
+  Eigen::RowVector3d p0(0., 0., 1.);
+
+  // Construct a sparse voxel grid whose cubes have edge length eps = 0.1.
+  // The cubes will form a thin shell around the implicit surface
+  const double eps = 0.1;
+
+  // CS will hold one scalar value at each cube vertex corresponding
+  // the value of the implicit at that vertex
+  Eigen::VectorXd CS;
+
+  // CV will hold the positions of the corners of the sparse voxel grid
+  Eigen::MatrixXd CV;
+
+  // CI is a #cubes x 8 matrix of indices where each row contains the
+  // indices into CV of the 8 corners of a cube
+  Eigen::MatrixXi CI;
+
+  // Construct the voxel grid, populating CS, CV, and CI
+  igl::sparse_voxel_grid(p0, scalar_func, eps, CS, CV, CI);
+
+  // Given the sparse voxel grid, use Marching Cubes to construct a triangle mesh of the surface
+  Eigen::MatrixXi F;
+  Eigen::MatrixXd V;
+  igl::copyleft::marching_cubes(CS, CV, CI, V, F);
+
+  // Draw the meshed implicit surface
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().clear();
+  viewer.data().set_mesh(V,F);
+  viewer.data().set_face_based(true);
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -142,6 +142,7 @@ if(TUTORIALS_CHAPTER7)
   if(LIBIGL_WITH_TETGEN)
     add_subdirectory("714_MarchingTets")
   endif()
+  add_subdirectory("715_MeshImplicitFunction")
 endif()