From 33d702375229874031a1418f052c70bd76cfeebd Mon Sep 17 00:00:00 2001 From: Marius Kintel Date: Sat, 14 Feb 2015 12:28:26 -0500 Subject: [PATCH] #1215 Repair non-manifold tessellation results, detect and report non-manifold results --- cgal/polyhole-tessellator-libtess2.cpp | 2 +- src/GeometryUtils.cc | 310 +++++++++++++++++++++---- src/GeometryUtils.h | 16 +- src/cgalutils.cc | 133 ++++++++++- 4 files changed, 409 insertions(+), 52 deletions(-) diff --git a/cgal/polyhole-tessellator-libtess2.cpp b/cgal/polyhole-tessellator-libtess2.cpp index f0ef286b..73bed56e 100644 --- a/cgal/polyhole-tessellator-libtess2.cpp +++ b/cgal/polyhole-tessellator-libtess2.cpp @@ -148,7 +148,7 @@ int main(int argc, char *argv[]) } std::vector triangles; - bool ok = GeometryUtils::tessellatePolygonWithHoles(polyhole, triangles, normal); + bool ok = GeometryUtils::tessellatePolygonWithHoles(&polyhole.vertices.front(), polyhole.faces, triangles, normal); std::cerr << "Tessellated into " << triangles.size() << " triangles" << std::endl; IndexedTriangleMesh trimesh; diff --git a/src/GeometryUtils.cc b/src/GeometryUtils.cc index 24a837e4..7f2191f2 100644 --- a/src/GeometryUtils.cc +++ b/src/GeometryUtils.cc @@ -4,6 +4,8 @@ #include "Reindexer.h" #include "grid.h" #include +#include +#include static void *stdAlloc(void* userData, unsigned int size) { TESS_NOTUSED(userData); @@ -15,6 +17,147 @@ static void stdFree(void* userData, void* ptr) { free(ptr); } +typedef std::pair IndexedEdge; + +/*! + Helper class for keeping track of edges in a mesh. + Can probably be replaced with a proper HalfEdge mesh later on +*/ +class EdgeDict { +public: +// Counts occurrences of edges + typedef boost::unordered_map IndexedEdgeDict; + + EdgeDict() { } + + void add(const IndexedFace &face) { + for (size_t i=0;icount(e) > 0) this->remove(e); + else this->add(e.second, e.first); + } + } + + void remove(const IndexedTriangle &t) { + for (int i=0;i<3;i++) { + IndexedEdge e(t[i], t[(i+1)%3]); + // If the edge exist, remove it + if (this->count(e) > 0) this->remove(e); + else this->add(e.second, e.first); + } + } + + void add(const IndexedTriangle &t) { + for (int i=0;i<3;i++) { + IndexedEdge e(t[(i+1)%3], t[i]); + // If an opposite edge exists, they cancel out + if (this->count(e) > 0) this->remove(e); + else this->add(e.second, e.first); + } + } + + void add(int start, int end) { + this->add(IndexedEdge(start,end)); + } + + void add(const IndexedEdge &e) { + this->edges[e]++; + PRINTDB("add: (%d,%d)", e.first % e.second); + } + + void remove(int start, int end) { + this->remove(IndexedEdge(start,end)); + } + + void remove(const IndexedEdge &e) { + this->edges[e]--; + if (this->edges[e] == 0) this->edges.erase(e); + PRINTDB("remove: (%d,%d)", e.first % e.second); + } + + int count(int start, int end) { + return this->count(IndexedEdge(start, end)); + } + + int count(const IndexedEdge &e) { + IndexedEdgeDict::const_iterator it = this->edges.find(e); + if (it != edges.end()) return it->second; + return 0; + } + + bool empty() const { return this->edges.empty(); } + + size_t size() const { return this->edges.size(); } + + void print() const { + BOOST_FOREACH(const IndexedEdgeDict::value_type &v, this->edges) { + const IndexedEdge &e = v.first; + PRINTDB(" (%d,%d)%s", e.first % e.second % ((v.second > 1) ? boost::lexical_cast(v.second).c_str() : "")); + } + } + + // Triangulate remaining loops and add to triangles + void triangulateLoops(std::vector &triangles) { + // First, look for self-intersections in edges + boost::unordered_map > v2e; + boost::unordered_map > v2e_reverse; + + BOOST_FOREACH(const IndexedEdgeDict::value_type &v, this->edges) { + const IndexedEdge &e = v.first; + v2e[e.first].push_back(e.second); + v2e_reverse[e.second].push_back(e.first); + } + + while (!v2e.empty()) { + for (boost::unordered_map >::iterator it = v2e.begin(); + it != v2e.end(); + it++) { + if (it->second.size() == 1) { // First single vertex + int vidx = it->first; + int next = it->second.front(); + assert(v2e_reverse.find(vidx) != v2e_reverse.end()); + assert(!v2e_reverse[vidx].empty()); + int prev = v2e_reverse[vidx].front(); + + IndexedTriangle t(prev, vidx, next); + PRINTDB("Clipping ear: %d %d %d", t[0] % t[1] % t[2]); + triangles.push_back(t); + // Remove the generated triangle from the original. + // Add new boundary edges to the edge dict + this->remove(t); + + v2e.erase(vidx); // single vertex + v2e_reverse.erase(vidx); // single vertex + + // If next->prev exists, remove it, otherwise add prev->next + std::list::iterator v2eit = std::find(v2e[next].begin(), v2e[next].end(), prev); + if (v2eit != v2e[next].end()) { + v2e[next].erase(v2eit); + v2e_reverse[prev].remove(next); + } + else { + v2e[prev].push_back(next); + v2e_reverse[next].push_back(prev); + } + if (v2e[next].empty()) v2e.erase(next); + if (v2e_reverse[prev].empty()) v2e.erase(prev); + + // Remove prev->vidx + v2e[prev].remove(vidx); + if (v2e[prev].empty()) v2e.erase(prev); + v2e_reverse[next].remove(vidx); + if (v2e_reverse[next].empty()) v2e_reverse.erase(next); + + break; + } + } + } + } + + IndexedEdgeDict edges; +}; + + /*! Tessellates input contours into a triangle mesh. @@ -33,23 +176,33 @@ static void stdFree(void* userData, void* ptr) { Returns true on error, false on success. */ -bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons, +bool GeometryUtils::tessellatePolygonWithHoles(const Vector3f *vertices, + const std::vector &faces, std::vector &triangles, const Vector3f *normal) { + // Algorithm outline: + // o Remove consecutive equal vertices and null ears (i.e. 23,24,23) + // o Ignore polygons with < 3 vertices + // o Pass-through polygons with 3 vertices + // o Pass polygon to libtess2 + // o Postprocess to clean up misbehaviors in libtess2 + // No polygon. FIXME: Will this ever happen or can we assert here? - if (polygons.faces.empty()) return false; + if (faces.empty()) return false; // Remove consecutive equal vertices, as well as null ears - std::vector faces = polygons.faces; - BOOST_FOREACH(IndexedFace &face, faces) { - int i=0; - while (i < face.size()) { + std::vector cleanfaces = faces; + BOOST_FOREACH(IndexedFace &face, cleanfaces) { + size_t i=0; + while (face.size() >= 3 && i < face.size()) { if (face[i] == face[(i+1)%face.size()]) { // Two consecutively equal indices face.erase(face.begin()+i); } - else if (face[i] == face[(i+2)%face.size()]) { // Null ear - face.erase(face.begin() + (i+1)%face.size()); + else if (face[(i+face.size()-1)%face.size()] == face[(i+1)%face.size()]) { // Null ear + if (i == 0) face.erase(face.begin() + i, face.begin() + i + 2); + else face.erase(face.begin() + i - 1, face.begin() + i + 1); + i--; } else { i++; @@ -57,22 +210,30 @@ bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons, } } // First polygon has < 3 points - no output - if (faces[0].size() < 3) return false; + if (cleanfaces[0].size() < 3) return false; // Remove collapsed holes - for (int i=1;i allindices; - BOOST_FOREACH(const IndexedFace &face, faces) { + BOOST_FOREACH(const IndexedFace &face, cleanfaces) { contour.clear(); BOOST_FOREACH(int idx, face) { - const Vector3f &v = verts[idx]; + const Vector3f &v = vertices[idx]; contour.push_back(v[0]); contour.push_back(v[1]); contour.push_back(v[2]); @@ -125,10 +286,10 @@ bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons, vertices for intersecting edges, we need to detect these and insert dummy triangles to maintain external connectivity. - FIXME: This currently only works for polygons without holes. + In addition, libtess2 may generate flipped triangles, i.e. triangles + where the edge direction is reversed compared to the input polygon. + This will also destroy connectivity and we need to flip those back. */ - if (faces.size() == 1) { // Only works for polygons without holes - /* Algorithm: A) Collect all triangles using _only_ existing vertices -> triangles @@ -162,17 +323,58 @@ bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons, vflags[tri[0]]++; // B) vflags[tri[1]]++; vflags[tri[2]]++; + + // For each edge in mappedtri, locate the opposite edge in the original polygon. + // If an opposite edge was found, we need to flip. + // In addition, remove each edge from the dict to be able to later find + // missing edges. + // Note: In some degenerate cases, we create triangles with mixed edge directions. + // In this case, don't reverse, but attempt to carry on + bool reverse = false; + for (int i=0;i<3;i++) { + const IndexedEdge e(mappedtri[i], mappedtri[(i+1)%3]); + if (edges.count(e) > 0) { + reverse = false; + break; + } + else if (edges.count(e.second, e.first) > 0) { + reverse = true; + } + } + if (reverse) { + mappedtri.reverseInPlace(); + PRINTDB(" reversed: %d %d %d", + mappedtri[0] % mappedtri[1] % mappedtri[2]); + } + + // Remove the generated triangle from the original. + // Add new boundary edges to the edge dict + edges.remove(mappedtri); triangles.push_back(mappedtri); } } + if (!edges.empty()) { + PRINTDB(" %d edges remaining after main triangulation", edges.size()); + edges.print(); + + // Collect loops from remaining edges and triangulate loops manually + edges.triangulateLoops(triangles); + + if (!edges.empty()) { + PRINTDB(" %d edges remaining after loop triangulation", edges.size()); + edges.print(); + } + } +#if 0 for (int i=0;i uniqueVertices; - IndexedPolygons indexedpolygons; - indexedpolygons.faces.push_back(IndexedFace()); - IndexedFace &currface = indexedpolygons.faces.back(); + std::vector indexedfaces; + indexedfaces.push_back(IndexedFace()); + IndexedFace &currface = indexedfaces.back(); BOOST_FOREACH (const Vector3d &v, polygon) { int idx = uniqueVertices.lookup(v.cast()); if (currface.empty() || idx != currface.back()) currface.push_back(idx); } if (currface.front() == currface.back()) currface.pop_back(); if (currface.size() >= 3) { // Cull empty triangles - uniqueVertices.copy(std::back_inserter(indexedpolygons.vertices)); + const Vector3f *verts = uniqueVertices.getArray(); std::vector indexedtriangles; - err = tessellatePolygonWithHoles(indexedpolygons, indexedtriangles, normal); - Vector3f *verts = &indexedpolygons.vertices.front(); + err = tessellatePolygonWithHoles(verts, indexedfaces, indexedtriangles, normal); BOOST_FOREACH(const IndexedTriangle &t, indexedtriangles) { triangles.push_back(Polygon()); Polygon &p = triangles.back(); @@ -242,3 +423,36 @@ bool GeometryUtils::tessellatePolygon(const Polygon &polygon, Polygons &triangle } return err; } + +int GeometryUtils::findUnconnectedEdges(const std::vector > &polygons) +{ + EdgeDict edges; + BOOST_FOREACH(const std::vector &faces, polygons) { + BOOST_FOREACH(const IndexedFace &face, faces) { + edges.add(face); + } + } +#if 1 // for debugging + if (!edges.empty()) { + PRINTD("Unconnected:"); + edges.print(); + } +#endif + return edges.size(); +} + +int GeometryUtils::findUnconnectedEdges(const std::vector &triangles) +{ + EdgeDict edges; + BOOST_FOREACH(const IndexedTriangle &t, triangles) { + edges.add(t); + } +#if 1 // for debugging + if (!edges.empty()) { + PRINTD("Unconnected:"); + edges.print(); + } +#endif + + return edges.size(); +} diff --git a/src/GeometryUtils.h b/src/GeometryUtils.h index 4d3b3f19..8d507adb 100644 --- a/src/GeometryUtils.h +++ b/src/GeometryUtils.h @@ -19,9 +19,21 @@ struct IndexedTriangleMesh { std::vector triangles; }; +// Indexed polygon mesh, where each polygon can have holes +struct IndexedPolyMesh { + std::vector vertices; + std::vector > polygons; +}; + namespace GeometryUtils { - bool tessellatePolygon(const Polygon &polygon, Polygons &triangles, + bool tessellatePolygon(const Polygon &polygon, + Polygons &triangles, const Vector3f *normal = NULL); - bool tessellatePolygonWithHoles(const IndexedPolygons &polygons, std::vector &triangles, + bool tessellatePolygonWithHoles(const Vector3f *vertices, + const std::vector &faces, + std::vector &triangles, const Vector3f *normal = NULL); + + int findUnconnectedEdges(const std::vector > &polygons); + int findUnconnectedEdges(const std::vector &triangles); } diff --git a/src/cgalutils.cc b/src/cgalutils.cc index a9466b54..6253d3de 100644 --- a/src/cgalutils.cc +++ b/src/cgalutils.cc @@ -1025,7 +1025,7 @@ namespace CGALUtils { return err; } #endif -#if 1 +#if 0 bool createPolySetFromNefPolyhedron3(const CGAL_Nef_polyhedron3 &N, PolySet &ps) { bool err = false; @@ -1092,6 +1092,137 @@ namespace CGALUtils { } } } + } +#endif +#if 1 + bool createPolySetFromNefPolyhedron3(const CGAL_Nef_polyhedron3 &N, PolySet &ps) + { + // 1. Build Indexed PolyMesh + // 2. Validate mesh (manifoldness) + // 3. Triangulate each face + // -> IndexedTriangleMesh + // 4. Validate mesh (manifoldness) + // 5. Create PolySet + + bool err = false; + + // 1. Build Indexed PolyMesh + Reindexer allVertices; + std::vector > polygons; + + CGAL_Nef_polyhedron3::Halffacet_const_iterator hfaceti; + CGAL_forall_halffacets(hfaceti, N) { + CGAL::Plane_3 plane(hfaceti->plane()); + // Since we're downscaling to float, vertices might merge during this conversion. + // To avoid passing equal vertices to the tessellator, we remove consecutively identical + // vertices. + polygons.push_back(std::vector()); + std::vector &faces = polygons.back(); + // the 0-mark-volume is the 'empty' volume of space. skip it. + if (!hfaceti->incident_volume()->mark()) { + CGAL_Nef_polyhedron3::Halffacet_cycle_const_iterator cyclei; + CGAL_forall_facet_cycles_of(cyclei, hfaceti) { + CGAL_Nef_polyhedron3::SHalfedge_around_facet_const_circulator c1(cyclei); + CGAL_Nef_polyhedron3::SHalfedge_around_facet_const_circulator c2(c1); + faces.push_back(IndexedFace()); + IndexedFace &currface = faces.back(); + CGAL_For_all(c1, c2) { + CGAL_Point_3 p = c1->source()->center_vertex()->point(); + // Create vertex indices and remove consecutive duplicate vertices + int idx = allVertices.lookup(vector_convert(p)); + if (currface.empty() || idx != currface.back()) currface.push_back(idx); + } + if (currface.front() == currface.back()) currface.pop_back(); + if (currface.size() < 3) faces.pop_back(); // Cull empty triangles + } + } + if (faces.empty()) polygons.pop_back(); // Cull empty faces + } + + // 2. Validate mesh (manifoldness) + int unconnected = GeometryUtils::findUnconnectedEdges(polygons); + if (unconnected > 0) { + PRINTB("Error: Non-manifold mesh encountered: %d unconnected edges", unconnected); + } + // 3. Triangulate each face + const Vector3f *verts = allVertices.getArray(); + std::vector allTriangles; + BOOST_FOREACH(const std::vector &faces, polygons) { +#if 0 // For debugging + std::cerr << "---\n"; + BOOST_FOREACH(const IndexedFace &poly, faces) { + BOOST_FOREACH(int i, poly) { + std::cerr << i << " "; + } + std::cerr << "\n"; + } +#if 0 + std::cerr.precision(20); + BOOST_FOREACH(const IndexedFace &poly, faces) { + BOOST_FOREACH(int i, poly) { + std::cerr << verts[i][0] << "," << verts[i][1] << "," << verts[i][2] << "\n"; + } + std::cerr << "\n"; + } +#endif + std::cerr << "-\n"; +#endif +#if 0 // For debugging + std::cerr.precision(20); + for (int i=0;i nvec = plane.orthogonal_vector(); + // K::Vector_3 normal(CGAL::to_double(nvec.x()), CGAL::to_double(nvec.y()), CGAL::to_double(nvec.z())); + std::vector triangles; + bool err = GeometryUtils::tessellatePolygonWithHoles(verts, faces, triangles, NULL); + if (!err) { + BOOST_FOREACH(const IndexedTriangle &t, triangles) { + assert(t[0] >= 0 && t[0] < allVertices.size()); + assert(t[1] >= 0 && t[1] < allVertices.size()); + assert(t[2] >= 0 && t[2] < allVertices.size()); + allTriangles.push_back(t); + } + } + } + +#if 0 // For debugging + BOOST_FOREACH(const IndexedTriangle &t, allTriangles) { + std::cerr << t[0] << " " << t[1] << " " << t[2] << "\n"; + } +#endif + // 4. Validate mesh (manifoldness) + int unconnected2 = GeometryUtils::findUnconnectedEdges(allTriangles); + if (unconnected2 > 0) { + PRINTB("Error: Non-manifold triangle mesh created: %d unconnected edges", unconnected2); + } + + BOOST_FOREACH(const IndexedTriangle &t, allTriangles) { + ps.append_poly(); + ps.append_vertex(verts[t[0]]); + ps.append_vertex(verts[t[1]]); + ps.append_vertex(verts[t[2]]); + } + +#if 0 // For debugging + std::cerr.precision(20); + for (int i=0;i