diff --git a/include/igl/copyleft/marching_cubes.h b/include/igl/copyleft/marching_cubes.h index 941d3dea2..925691a67 100644 --- a/include/igl/copyleft/marching_cubes.h +++ b/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. // diff --git a/include/igl/sparse_voxel_grid.cpp b/include/igl/sparse_voxel_grid.cpp new file mode 100644 index 000000000..1ffb8bc0b --- /dev/null +++ b/include/igl/sparse_voxel_grid.cpp @@ -0,0 +1,140 @@ +#include "sparse_voxel_grid.h" + +#include +#include +#include + + +template +IGL_INLINE void igl::sparse_voxel_grid(const Eigen::MatrixBase& p0, + const std::function& scalarFunc, + const double eps, + Eigen::PlainObjectBase& CS, + Eigen::PlainObjectBase& CV, + Eigen::PlainObjectBase& CI, + int expected_number_of_cubes) +{ + typedef typename DerivedV::Scalar ScalarV; + typedef typename DerivedS::Scalar ScalarS; + typedef typename DerivedI::Scalar ScalarI; + typedef Eigen::Matrix VertexRowVector; + typedef Eigen::Matrix IndexRowVector; + + + struct IndexRowVectorHash { + std::size_t operator()(const Eigen::RowVector3i& key) const { + std::size_t seed = 0; + std::hash 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 CI_vector(expected_number_of_cubes); + std::vector CV_vector(8*expected_number_of_cubes); + std::vector CS_vector(8*expected_number_of_cubes); + + // Track visisted neighbors + std::unordered_map visited(6*expected_number_of_cubes); + + // BFS Queue + std::vector 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(); // 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 neighbors = { + bx, -bx, by, -by, bz, -bz + }; + + // Compute the position of the cube corners and the scalar values at those corners + std::array cubeCorners = { + ctr+half_eps*(bx+by+bz).cast(), ctr+half_eps*(bx+by-bz).cast(), ctr+half_eps*(-bx+by-bz).cast(), ctr+half_eps*(-bx+by+bz).cast(), + ctr+half_eps*(bx-by+bz).cast(), ctr+half_eps*(bx-by-bz).cast(), ctr+half_eps*(-bx-by-bz).cast(), ctr+half_eps*(-bx-by+bz).cast() + }; + std::array 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 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, 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, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, std::function::Scalar (Eigen::Matrix const&)> const&, double, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, int); +#endif diff --git a/include/igl/sparse_voxel_grid.h b/include/igl/sparse_voxel_grid.h new file mode 100644 index 000000000..da1fcd940 --- /dev/null +++ b/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 +// +// 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 + +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 + IGL_INLINE void sparse_voxel_grid(const Eigen::MatrixBase& p0, + const std::function& scalarFunc, + const double eps, + Eigen::PlainObjectBase& CS, + Eigen::PlainObjectBase& CV, + Eigen::PlainObjectBase& CI, + int expected_number_of_cubes=1024); + +} +#ifndef IGL_STATIC_LIBRARY +# include "sparse_voxel_grid.cpp" +#endif + +#endif // IGL_SPARSE_VOXEL_GRID_H diff --git a/tutorial/715_MeshImplicitFunction/CMakeLists.txt b/tutorial/715_MeshImplicitFunction/CMakeLists.txt new file mode 100644 index 000000000..c4ce2b52d --- /dev/null +++ b/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) diff --git a/tutorial/715_MeshImplicitFunction/main.cpp b/tutorial/715_MeshImplicitFunction/main.cpp new file mode 100755 index 000000000..ae693151c --- /dev/null +++ b/tutorial/715_MeshImplicitFunction/main.cpp @@ -0,0 +1,50 @@ +#include +#include +#include + +#include +#include + +#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 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(); +} diff --git a/tutorial/CMakeLists.txt b/tutorial/CMakeLists.txt index 94dd96220..5f8b7a6f6 100644 --- a/tutorial/CMakeLists.txt +++ b/tutorial/CMakeLists.txt @@ -142,6 +142,7 @@ if(TUTORIALS_CHAPTER7) if(LIBIGL_WITH_TETGEN) add_subdirectory("714_MarchingTets") endif() + add_subdirectory("715_MeshImplicitFunction") endif()