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.
This commit is contained in:
Francis Williams
committed by
Daniele Panozzo
Daniele Panozzo
parent
fa773a79f4
commit
4f030acb67
@@ -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.
|
||||
//
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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)
|
||||
Executable
+50
@@ -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();
|
||||
}
|
||||
@@ -142,6 +142,7 @@ if(TUTORIALS_CHAPTER7)
|
||||
if(LIBIGL_WITH_TETGEN)
|
||||
add_subdirectory("714_MarchingTets")
|
||||
endif()
|
||||
add_subdirectory("715_MeshImplicitFunction")
|
||||
endif()
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user