diff --git a/include/igl/adjacency_list.cpp b/include/igl/adjacency_list.cpp index 1fb61ebf1..dfc94be2f 100644 --- a/include/igl/adjacency_list.cpp +++ b/include/igl/adjacency_list.cpp @@ -165,4 +165,5 @@ template void igl::adjacency_list, int>(Eige // generated by autoexplicit.sh template void igl::adjacency_list, int>(Eigen::PlainObjectBase > const&, std::vector >, std::allocator > > >&, bool); template void igl::adjacency_list, int>(Eigen::PlainObjectBase > const&, std::vector >, std::allocator > > >&, bool); +template void igl::adjacency_list, unsigned int>(class Eigen::PlainObjectBase > const &, class std::vector >, class std::allocator > > > &, bool); #endif diff --git a/include/igl/centroid.cpp b/include/igl/centroid.cpp index a9ed25df2..7587aeddd 100644 --- a/include/igl/centroid.cpp +++ b/include/igl/centroid.cpp @@ -64,4 +64,5 @@ IGL_INLINE void igl::centroid( // generated by autoexplicit.sh template void igl::centroid, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&); template void igl::centroid, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&); +template void igl::centroid, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&); #endif diff --git a/include/igl/copyleft/marching_cubes.cpp b/include/igl/copyleft/marching_cubes.cpp index 144a02255..c385cbb70 100644 --- a/include/igl/copyleft/marching_cubes.cpp +++ b/include/igl/copyleft/marching_cubes.cpp @@ -33,44 +33,45 @@ extern const int edgeTable[256]; extern const int triTable[256][2][17]; extern const int polyTable[8][16]; -struct EdgeKey -{ - EdgeKey(unsigned i0, unsigned i1) : i0_(i0), i1_(i1) {} - bool operator==(const EdgeKey& _rhs) const - { - return i0_ == _rhs.i0_ && i1_ == _rhs.i1_; - } - - unsigned i0_, i1_; -}; - -struct EdgeHash -{ - std::size_t operator()(const EdgeKey& key) const { - std::size_t seed = 0; - seed ^= key.i0_ + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine - seed ^= key.i1_ + 0x9e3779b9 + (seed<<6) + (seed>>2); - return std::hash()(seed); - } -}; - - -template +template class MarchingCubes { + struct EdgeKey + { + EdgeKey(unsigned i0, unsigned i1) : i0_(i0), i1_(i1) {} + + bool operator==(const EdgeKey& _rhs) const + { + return i0_ == _rhs.i0_ && i1_ == _rhs.i1_; + } + + unsigned i0_, i1_; + }; + + struct EdgeHash + { + std::size_t operator()(const EdgeKey& key) const { + std::size_t seed = 0; + seed ^= key.i0_ + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine + seed ^= key.i1_ + 0x9e3779b9 + (seed<<6) + (seed>>2); + return std::hash()(seed); + } + }; + typedef std::unordered_map MyMap; typedef typename MyMap::const_iterator MyMapIterator; public: - MarchingCubes( - const Eigen::PlainObjectBase &values, - const Eigen::PlainObjectBase &points, + // Dense index grid version + MarchingCubes(const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, const unsigned x_res, const unsigned y_res, const unsigned z_res, - Eigen::PlainObjectBase &vertices, - Eigen::PlainObjectBase &faces) + const double isovalue, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces) { assert(values.cols() == 1); assert(points.cols() == 3); @@ -101,7 +102,7 @@ public: { unsigned corner[8]; - typename DerivedF::Scalar samples[12]; + typename DerivedFaces::Scalar samples[12]; unsigned char cubetype(0); unsigned int i; @@ -126,7 +127,7 @@ public: // determine cube type for (i=0; i<8; ++i) - if (values(corner[i]) > 0.0) + if (values(corner[i]) > isovalue) cubetype |= (1< &values, - const Eigen::PlainObjectBase &points, - unsigned int i0, - unsigned int i1, - Eigen::PlainObjectBase &vertices, - int &num_vertices, - MyMap &edge2vertex) + // Sparse index grid version + MarchingCubes(const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const Eigen::MatrixBase &cubes, + const double isovalue, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces) + { + assert(values.cols() == 1); + assert(points.cols() == 3); + assert(cubes.cols() == 8); + + if(cubes.rows() == 0) + { + return; + } + + faces.resize(10000,3); + int num_faces = 0; + + vertices.resize(10000,3); + int num_vertices = 0; + + + unsigned n_cubes = cubes.rows(); + + for (unsigned cube_it =0 ; cube_it < n_cubes; ++cube_it) + { + typedef Eigen::Matrix CubeIndexVector; + typedef typename DerivedFaces::Scalar SampleScalar; + + CubeIndexVector cube = cubes.row(cube_it); + SampleScalar samples[12]; + unsigned char cubetype(0); + + // determine cube type + for (int i=0; i<8; ++i) + { + if (values[cube[i]] > isovalue) + { + cubetype |= (1< faces.rows()) + { + faces.conservativeResize(faces.rows()+10000, Eigen::NoChange); + } + faces.row(num_faces-1) << + samples[triTable[cubetype][0][i ]], + samples[triTable[cubetype][0][i+1]], + samples[triTable[cubetype][0][i+2]]; + + } + + } + + vertices.conservativeResize(num_vertices, Eigen::NoChange); + faces.conservativeResize(num_faces, Eigen::NoChange); + + } + + static typename DerivedFaces::Scalar add_vertex(const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + unsigned int i0, + unsigned int i1, + Eigen::PlainObjectBase &vertices, + int &num_vertices, + MyMap &edge2vertex) { // find vertex if it has been computed already MyMapIterator it = edge2vertex.find(EdgeKey(i0, i1)); if (it != edge2vertex.end()) + { return it->second; - ; + } // generate new vertex - const Eigen::Matrix & p0 = points.row(i0); - const Eigen::Matrix & p1 = points.row(i1); - - typename Derivedvalues::Scalar s0 = fabs(values(i0)); - typename Derivedvalues::Scalar s1 = fabs(values(i1)); - typename Derivedvalues::Scalar t = s0 / (s0+s1); + const Eigen::Matrix & p0 = points.row(i0); + const Eigen::Matrix & p1 = points.row(i1); + typename DerivedValues::Scalar s0 = fabs(values[i0]); + typename DerivedValues::Scalar s1 = fabs(values[i1]); + typename DerivedValues::Scalar t = s0 / (s0+s1); num_vertices++; if (num_vertices > vertices.rows()) + { vertices.conservativeResize(vertices.rows()+10000, Eigen::NoChange); + } // Linear interpolation based on linearly interpolating values - vertices.row(num_vertices-1) = ((1.0f-t)*p0 + t*p1).template cast(); + vertices.row(num_vertices-1) = ((1.0f-t)*p0 + t*p1).template cast(); edge2vertex[EdgeKey(i0, i1)] = num_vertices-1; return num_vertices-1; } - ; // maps an edge to the sample vertex generated on it MyMap edge2vertex; }; -template +template IGL_INLINE void igl::copyleft::marching_cubes( - const Eigen::PlainObjectBase &values, - const Eigen::PlainObjectBase &points, + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, const unsigned x_res, const unsigned y_res, const unsigned z_res, - Eigen::PlainObjectBase &vertices, - Eigen::PlainObjectBase &faces) + const double isovalue, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces) { - MarchingCubes mc(values, - points, - x_res, - y_res, - z_res, - vertices, - faces); + typedef Eigen::MatrixXi Shim; /* DerivedIndices shim type is unused in this instantiation*/ + MarchingCubes + mc(values, points, x_res, y_res, z_res, isovalue, vertices, faces); } + +template +IGL_INLINE void igl::copyleft::marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const unsigned x_res, + const unsigned y_res, + const unsigned z_res, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces) +{ + typedef Eigen::MatrixXi Shim; /* DerivedIndices shim type is unused in this instantiation*/ + MarchingCubes + mc(values, points, x_res, y_res, z_res, 0.0 /*isovalue*/, vertices, faces); +} + +template +IGL_INLINE void igl::copyleft::marching_cubes( + const Eigen::MatrixBase& values, + const Eigen::MatrixBase& points, + const Eigen::MatrixBase& indices, + const double isovalue, + Eigen::PlainObjectBase& vertices, + Eigen::PlainObjectBase &faces) +{ + MarchingCubes mc(values, points, indices, isovalue, vertices, faces); +} + +template +IGL_INLINE void igl::copyleft::marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const Eigen::MatrixBase & indices, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces) +{ + MarchingCubes mc(values, points, indices, 0.0 /*isovalue*/, vertices, faces); +} + #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation + // generated by autoexplicit.sh -template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); // generated by autoexplicit.sh -template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); // generated by autoexplicit.sh -template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); // generated by autoexplicit.sh -template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); -// generated by autoexplicit.sh -template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); -template void igl::copyleft::marching_cubes< Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, unsigned int, unsigned int, unsigned int, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, const Eigen::MatrixBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix >(const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,Eigen::PlainObjectBase >&,Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix >(const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,Eigen::PlainObjectBase >&,Eigen::PlainObjectBase >&); +template void igl::copyleft::marching_cubes,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix,Eigen::Matrix >(const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,const Eigen::MatrixBase >&,Eigen::PlainObjectBase >&,Eigen::PlainObjectBase >&); #endif diff --git a/include/igl/copyleft/marching_cubes.h b/include/igl/copyleft/marching_cubes.h index 4c9f888e4..941d3dea2 100644 --- a/include/igl/copyleft/marching_cubes.h +++ b/include/igl/copyleft/marching_cubes.h @@ -14,10 +14,10 @@ namespace igl { namespace copyleft { - // marching_cubes( values, points, x_res, y_res, z_res, vertices, faces ) + // marching_cubes( values, points, x_res, y_res, z_res, isovalue, vertices, faces ) // - // performs marching cubes reconstruction on the grid defined by values, and - // points, and generates vertices and faces + // performs marching cubes reconstruction on a grid defined by values, and + // points, and generates a mesh defined by vertices and faces // // Input: // values #number_of_grid_points x 1 array -- the scalar values of an @@ -34,24 +34,75 @@ namespace igl // xres resolutions of the grid in x dimension // yres resolutions of the grid in y dimension // zres resolutions of the grid in z dimension + // isovalue the isovalue of the surface to reconstruct // Output: // vertices #V by 3 list of mesh vertex positions // faces #F by 3 list of mesh triangle indices // - template < - typename Derivedvalues, - typename Derivedpoints, - typename Derivedvertices, - typename DerivedF> - IGL_INLINE void marching_cubes( - const Eigen::PlainObjectBase &values, - const Eigen::PlainObjectBase &points, + template + IGL_INLINE void marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, const unsigned x_res, const unsigned y_res, const unsigned z_res, - Eigen::PlainObjectBase &vertices, - Eigen::PlainObjectBase &faces); + const double isovalue, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces); + + // Overload of the above function where the isovalue defaults to 0.0 + template + IGL_INLINE void marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const unsigned x_res, + const unsigned y_res, + const unsigned z_res, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces); + + + + // marching_cubes( values, points, indices, vertices, faces ) + // + // Perform marching cubes reconstruction on the grid cells defined by indices. + // 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. + // + // Input: + // values #number_of_grid_points x 1 array -- the scalar values of an + // implicit function defined on the grid points (<0 in the inside of the + // surface, 0 on the border, >0 outside) + // points #number_of_grid_points x 3 array -- 3-D positions of the grid + // points, ordered in x,y,z order: + // indices #cubes x 8 array -- one row for each cube where each value is + // the index of a vertex in points and a scalar in values. + // i.e. points[indices[i, j]] = the position of the j'th vertex of the i'th cube + // Output: + // vertices #V by 3 list of mesh vertex positions + // faces #F by 3 list of mesh triangle indices + // Note: The winding direction of the cube indices will affect the output winding of the faces + // + template + IGL_INLINE void marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const Eigen::MatrixBase &indices, + const double isovalue, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces); + + // Overload of the above function where the isovalue defaults to 0.0 + template + IGL_INLINE void marching_cubes( + const Eigen::MatrixBase &values, + const Eigen::MatrixBase &points, + const Eigen::MatrixBase &indices, + Eigen::PlainObjectBase &vertices, + Eigen::PlainObjectBase &faces); + } + } #ifndef IGL_STATIC_LIBRARY diff --git a/include/igl/cut_to_disk.cpp b/include/igl/cut_to_disk.cpp new file mode 100644 index 000000000..931c26714 --- /dev/null +++ b/include/igl/cut_to_disk.cpp @@ -0,0 +1,330 @@ +#include "cut_to_disk.h" + +#include +#include +#include +#include + +namespace igl { + template + void cut_to_disk( + const Eigen::MatrixBase &F, + std::vector > &cuts) + { + cuts.clear(); + + Index nfaces = F.rows(); + + if (nfaces == 0) + return; + + std::map, std::vector > edges; + // build edges + + for (Index i = 0; i < nfaces; i++) + { + for (int j = 0; j < 3; j++) + { + Index v0 = F(i, j); + Index v1 = F(i, (j + 1) % 3); + std::pair e; + e.first = std::min(v0, v1); + e.second = std::max(v0, v1); + edges[e].push_back(i); + } + } + + int nedges = edges.size(); + Eigen::Matrix edgeVerts(nedges,2); + Eigen::Matrix edgeFaces(nedges,2); + Eigen::Matrix faceEdges(nfaces, 3); + std::set boundaryEdges; + std::map, Index> edgeidx; + Index idx = 0; + for (auto it : edges) + { + edgeidx[it.first] = idx; + edgeVerts(idx, 0) = it.first.first; + edgeVerts(idx, 1) = it.first.second; + edgeFaces(idx, 0) = it.second[0]; + if (it.second.size() > 1) + { + edgeFaces(idx, 1) = it.second[1]; + } + else + { + edgeFaces(idx, 1) = -1; + boundaryEdges.insert(idx); + } + idx++; + } + for (Index i = 0; i < nfaces; i++) + { + for (int j = 0; j < 3; j++) + { + Index v0 = F(i, j); + Index v1 = F(i, (j + 1) % 3); + std::pair e; + e.first = std::min(v0, v1); + e.second = std::max(v0, v1); + faceEdges(i, j) = edgeidx[e]; + } + } + + bool *deleted = new bool[nfaces]; + for (Index i = 0; i < nfaces; i++) + deleted[i] = false; + + std::set deletededges; + + // loop over faces + for (Index face = 0; face < nfaces; face++) + { + // stop at first undeleted face + if (deleted[face]) + continue; + deleted[face] = true; + std::deque processEdges; + for (int i = 0; i < 3; i++) + { + Index e = faceEdges(face, i); + if (boundaryEdges.count(e)) + continue; + int ndeleted = 0; + if (deleted[edgeFaces(e, 0)]) + ndeleted++; + if (deleted[edgeFaces(e, 1)]) + ndeleted++; + if (ndeleted == 1) + processEdges.push_back(e); + } + // delete all faces adjacent to edges with exactly one adjacent face + while (!processEdges.empty()) + { + Index nexte = processEdges.front(); + processEdges.pop_front(); + Index todelete = nfaces; + if (!deleted[edgeFaces(nexte, 0)]) + todelete = edgeFaces(nexte, 0); + if (!deleted[edgeFaces(nexte, 1)]) + todelete = edgeFaces(nexte, 1); + if (todelete != nfaces) + { + deletededges.insert(nexte); + deleted[todelete] = true; + for (int i = 0; i < 3; i++) + { + Index e = faceEdges(todelete, i); + if (boundaryEdges.count(e)) + continue; + int ndeleted = 0; + if (deleted[edgeFaces(e, 0)]) + ndeleted++; + if (deleted[edgeFaces(e, 1)]) + ndeleted++; + if (ndeleted == 1) + processEdges.push_back(e); + } + } + } + } + delete[] deleted; + + // accumulated non-deleted edges + std::vector leftedges; + for (Index i = 0; i < nedges; i++) + { + if (!deletededges.count(i)) + leftedges.push_back(i); + } + + deletededges.clear(); + // prune spines + std::map > spinevertedges; + for (Index i : leftedges) + { + spinevertedges[edgeVerts(i, 0)].push_back(i); + spinevertedges[edgeVerts(i, 1)].push_back(i); + } + + std::deque vertsProcess; + std::map spinevertnbs; + for (auto it : spinevertedges) + { + spinevertnbs[it.first] = it.second.size(); + if (it.second.size() == 1) + vertsProcess.push_back(it.first); + } + while (!vertsProcess.empty()) + { + Index vert = vertsProcess.front(); + vertsProcess.pop_front(); + for (Index e : spinevertedges[vert]) + { + if (!deletededges.count(e)) + { + deletededges.insert(e); + for (int j = 0; j < 2; j++) + { + spinevertnbs[edgeVerts(e, j)]--; + if (spinevertnbs[edgeVerts(e, j)] == 1) + { + vertsProcess.push_back(edgeVerts(e, j)); + } + } + } + } + } + std::vector loopedges; + for (Index i : leftedges) + if (!deletededges.count(i)) + loopedges.push_back(i); + + Index nloopedges = loopedges.size(); + if (nloopedges == 0) + return; + + std::map > loopvertedges; + for (Index e : loopedges) + { + loopvertedges[edgeVerts(e, 0)].push_back(e); + loopvertedges[edgeVerts(e, 1)].push_back(e); + } + + std::set usededges; + for (Index e : loopedges) + { + // make a cycle or chain starting from this edge + while (!usededges.count(e)) + { + std::vector cycleverts; + std::vector cycleedges; + cycleverts.push_back(edgeVerts(e, 0)); + cycleverts.push_back(edgeVerts(e, 1)); + cycleedges.push_back(e); + + std::map cycleidx; + cycleidx[cycleverts[0]] = 0; + cycleidx[cycleverts[1]] = 1; + + Index curvert = edgeVerts(e, 1); + Index cure = e; + bool foundcycle = false; + while (curvert != -1 && !foundcycle) + { + Index nextvert = -1; + Index nexte = -1; + for (Index cande : loopvertedges[curvert]) + { + if (!usededges.count(cande) && cande != cure) + { + int vidx = 0; + if (curvert == edgeVerts(cande, vidx)) + vidx = 1; + nextvert = edgeVerts(cande, vidx); + nexte = cande; + break; + } + } + if (nextvert != -1) + { + auto it = cycleidx.find(nextvert); + if (it != cycleidx.end()) + { + // we've hit outselves + std::vector cut; + for (Index i = it->second; i < cycleverts.size(); i++) + { + cut.push_back(cycleverts[i]); + } + cut.push_back(nextvert); + cuts.push_back(cut); + for (Index i = it->second; i < cycleedges.size(); i++) + { + usededges.insert(cycleedges[i]); + } + usededges.insert(nexte); + foundcycle = true; + } + else + { + cycleidx[nextvert] = cycleverts.size(); + cycleverts.push_back(nextvert); + cycleedges.push_back(nexte); + } + } + curvert = nextvert; + cure = nexte; + } + if (!foundcycle) + { + // we've hit a dead end. reverse and try the other direction + std::reverse(cycleverts.begin(), cycleverts.end()); + std::reverse(cycleedges.begin(), cycleedges.end()); + curvert = cycleverts.back(); + cure = cycleedges.back(); + while (curvert != -1 && !foundcycle) + { + Index nextvert = -1; + Index nexte = -1; + for (Index cande : loopvertedges[curvert]) + { + if (!usededges.count(cande) && cande != cure) + { + int vidx = 0; + if (curvert == edgeVerts(cande, vidx)) + vidx = 1; + nextvert = edgeVerts(cande, vidx); + nexte = cande; + break; + } + } + if (nextvert != -1) + { + auto it = cycleidx.find(nextvert); + if (it != cycleidx.end()) + { + // we've hit outselves + std::vector cut; + for (Index i = it->second; i < cycleverts.size(); i++) + { + cut.push_back(cycleverts[i]); + } + cut.push_back(nextvert); + cuts.push_back(cut); + for (Index i = it->second; i < cycleedges.size(); i++) + { + usededges.insert(cycleedges[i]); + } + usededges.insert(nexte); + foundcycle = true; + } + else + { + cycleidx[nextvert] = cycleverts.size(); + cycleverts.push_back(nextvert); + cycleedges.push_back(nexte); + } + } + curvert = nextvert; + cure = nexte; + } + if (!foundcycle) + { + // we've found a chain + std::vector cut; + for (Index i = 0; i < cycleverts.size(); i++) + { + cut.push_back(cycleverts[i]); + } + cuts.push_back(cut); + for (Index i = 0; i < cycleedges.size(); i++) + { + usededges.insert(cycleedges[i]); + } + } + } + } + } + } +} diff --git a/include/igl/cut_to_disk.h b/include/igl/cut_to_disk.h new file mode 100644 index 000000000..cc6ae3266 --- /dev/null +++ b/include/igl/cut_to_disk.h @@ -0,0 +1,52 @@ +#ifndef IGL_CUT_TO_DISK_H +#define IGL_CUT_TO_DISK_H +#include "igl_inline.h" + +#include + +#include + +namespace igl +{ + // Given a triangle mesh, computes a set of edge cuts sufficient to carve the + // mesh into a topological disk, without disconnecting any connected components. + // Nothing else about the cuts (including number, total length, or smoothness) + // is guaranteed to be optimal. + // + // Simply-connected components without boundary (topological spheres) are left + // untouched (delete any edge if you really want a disk). + // All other connected components are cut into disks. Meshes with boundary are + // supported; boundary edges will be included as cuts. + // + // The cut mesh itself can be materialized using cut_mesh(). + // + // Implements the triangle-deletion approach described by Gu et al's + // "Geometry Images." + // + // Template Parameters: + // Index Integrable type large enough to represent the total number of faces + // and edges in the surface represented by F, and all entries of F. + // + // Inputs: + // F #F by 3 list of the faces (must be triangles) + // + // Outputs: + // cuts List of cuts. Each cut is a sequence of vertex indices (where + // pairs of consecutive vertices share a face), is simple, and is either + // a closed loop (in which the first and last indices are identical) or + // an open curve. Cuts are edge-disjoint. + // + + template < + typename DerivedF, + typename Index> + IGL_INLINE void cut_to_disk( + const Eigen::MatrixBase &F, + std::vector > &cuts); +}; + +#ifndef IGL_STATIC_LIBRARY +#include "cut_to_disk.cpp" +#endif + +#endif diff --git a/include/igl/dijkstra.cpp b/include/igl/dijkstra.cpp index ad92f1d7f..5c49936a1 100644 --- a/include/igl/dijkstra.cpp +++ b/include/igl/dijkstra.cpp @@ -5,14 +5,16 @@ // 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/. -#include +#include "dijkstra.h" template -IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source, - const std::set &targets, - const std::vector >& VV, - Eigen::PlainObjectBase &min_distance, - Eigen::PlainObjectBase &previous) +IGL_INLINE int igl::dijkstra( + const IndexType &source, + const std::set &targets, + const std::vector >& VV, + const std::vector& weights, + Eigen::PlainObjectBase &min_distance, + Eigen::PlainObjectBase &previous) { int numV = VV.size(); min_distance.setConstant(numV, 1, std::numeric_limits::infinity()); @@ -37,7 +39,7 @@ IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source, neighbor_iter++) { IndexType v = *neighbor_iter; - typename DerivedD::Scalar distance_through_u = dist + 1.; + typename DerivedD::Scalar distance_through_u = dist + weights[u]; if (distance_through_u < min_distance[v]) { vertex_queue.erase(std::make_pair(min_distance[v], v)); @@ -53,10 +55,23 @@ IGL_INLINE int igl::dijkstra_compute_paths(const IndexType &source, return -1; } +template +IGL_INLINE int igl::dijkstra( + const IndexType &source, + const std::set &targets, + const std::vector >& VV, + Eigen::PlainObjectBase &min_distance, + Eigen::PlainObjectBase &previous) +{ + std::vector weights(VV.size(), 1.0); + return dijkstra(source, targets, VV, weights, min_distance, previous); +} + template -IGL_INLINE void igl::dijkstra_get_shortest_path_to(const IndexType &vertex, - const Eigen::PlainObjectBase &previous, - std::vector &path) +IGL_INLINE void igl::dijkstra( + const IndexType &vertex, + const Eigen::MatrixBase &previous, + std::vector &path) { IndexType source = vertex; path.clear(); @@ -66,6 +81,6 @@ IGL_INLINE void igl::dijkstra_get_shortest_path_to(const IndexType &vertex, #ifdef IGL_STATIC_LIBRARY // Explicit template instantiation -template int igl::dijkstra_compute_paths, Eigen::Matrix >(int const&, std::set, std::allocator > const&, std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); -template void igl::dijkstra_get_shortest_path_to >(int const&, Eigen::PlainObjectBase > const&, std::vector >&); +template int igl::dijkstra, Eigen::Matrix >(int const&, std::set, std::allocator > const&, std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::dijkstra >(int const&, Eigen::MatrixBase > const&, std::vector >&); #endif diff --git a/include/igl/dijkstra.h b/include/igl/dijkstra.h index 33af160bb..933b426ce 100644 --- a/include/igl/dijkstra.h +++ b/include/igl/dijkstra.h @@ -16,7 +16,30 @@ namespace igl { - // Dijstra's algorithm for shortest paths, with multiple targets. + // Dijkstra's algorithm for vertex-weighted shortest paths, with multiple targets. + // Adapted from http://rosettacode.org/wiki/Dijkstra%27s_algorithm . + // + // Inputs: + // source index of source vertex + // targets target vector set + // VV #V list of lists of incident vertices (adjacency list), e.g. + // as returned by igl::adjacency_list + // weights #V list of scalar vertex weights + // + // Output: + // min_distance #V by 1 list of the minimum distances from source to all vertices + // previous #V by 1 list of the previous visited vertices (for each vertex) - used for backtracking + // + template + IGL_INLINE int dijkstra( + const IndexType &source, + const std::set &targets, + const std::vector >& VV, + const std::vector& weights, + Eigen::PlainObjectBase &min_distance, + Eigen::PlainObjectBase &previous); + + // Dijkstra's algorithm for shortest paths, with multiple targets. // Adapted from http://rosettacode.org/wiki/Dijkstra%27s_algorithm . // // Inputs: @@ -30,13 +53,14 @@ namespace igl { // previous #V by 1 list of the previous visited vertices (for each vertex) - used for backtracking // template - IGL_INLINE int dijkstra_compute_paths(const IndexType &source, - const std::set &targets, - const std::vector >& VV, - Eigen::PlainObjectBase &min_distance, - Eigen::PlainObjectBase &previous); + IGL_INLINE int dijkstra( + const IndexType &source, + const std::set &targets, + const std::vector >& VV, + Eigen::PlainObjectBase &min_distance, + Eigen::PlainObjectBase &previous); - // Backtracking after Dijstra's algorithm, to find shortest path. + // Backtracking after Dijkstra's algorithm, to find shortest path. // // Inputs: // vertex vertex to which we want the shortest path (from same source as above) @@ -46,9 +70,10 @@ namespace igl { // path #P by 1 list of vertex indices in the shortest path from source to vertex // template - IGL_INLINE void dijkstra_get_shortest_path_to(const IndexType &vertex, - const Eigen::PlainObjectBase &previous, - std::vector &path); + IGL_INLINE void dijkstra( + const IndexType &vertex, + const Eigen::MatrixBase &previous, + std::vector &path); }; diff --git a/include/igl/knn.cpp b/include/igl/knn.cpp index b2e355df1..5a415d3fa 100644 --- a/include/igl/knn.cpp +++ b/include/igl/knn.cpp @@ -102,4 +102,7 @@ namespace igl { +#ifdef IGL_STATIC_LIBRARY +template void igl::knn, int, int, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, int const&, std::vector >, std::allocator > > > const&, Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&); +#endif diff --git a/include/igl/polygon_mesh_to_triangle_mesh.cpp b/include/igl/polygon_mesh_to_triangle_mesh.cpp index 2ea43238b..ef10dc9dc 100644 --- a/include/igl/polygon_mesh_to_triangle_mesh.cpp +++ b/include/igl/polygon_mesh_to_triangle_mesh.cpp @@ -73,4 +73,5 @@ template void igl::polygon_mesh_to_triangle_mesh >(std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&); template void igl::polygon_mesh_to_triangle_mesh >(std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&); template void igl::polygon_mesh_to_triangle_mesh >(std::vector >, std::allocator > > > const&, Eigen::PlainObjectBase >&); +template void igl::polygon_mesh_to_triangle_mesh, Eigen::Matrix >(Eigen::PlainObjectBase > const&, Eigen::PlainObjectBase >&); #endif diff --git a/include/igl/remove_unreferenced.cpp b/include/igl/remove_unreferenced.cpp index 3ee8260aa..aea9f48ba 100644 --- a/include/igl/remove_unreferenced.cpp +++ b/include/igl/remove_unreferenced.cpp @@ -119,4 +119,5 @@ template void igl::remove_unreferenced, E template void igl::remove_unreferenced, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::remove_unreferenced, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); template void igl::remove_unreferenced, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); +template void igl::remove_unreferenced, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&, Eigen::PlainObjectBase >&); #endif