#1215 Repair non-manifold tessellation results, detect and report non-manifold results

cgal/polyhole-tessellator-libtess2.cpp

@@ -148,7 +148,7 @@ int main(int argc, char *argv[])
   }
 
   std::vector<IndexedTriangle> 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;

src/GeometryUtils.cc

@@ -4,6 +4,8 @@
 #include "Reindexer.h"
 #include "grid.h"
 #include <boost/foreach.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/unordered_map.hpp>
 
 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<int,int> 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<IndexedEdge, int> IndexedEdgeDict;
+
+	EdgeDict() { }
+
+	void add(const IndexedFace &face) {
+		for (size_t i=0;i<face.size();i++) {
+			IndexedEdge e(face[(i+1)%face.size()], face[i]);
+			if (this->count(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<std::string>(v.second).c_str() : ""));
+		}
+	}
+
+	// Triangulate remaining loops and add to triangles
+	void triangulateLoops(std::vector<IndexedTriangle> &triangles) {
+		// First, look for self-intersections in edges
+		boost::unordered_map<int, std::list<int> > v2e;
+		boost::unordered_map<int, std::list<int> > 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<int, std::list<int> >::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<int>::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<IndexedFace> &faces, 
 																							 std::vector<IndexedTriangle> &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<IndexedFace> faces = polygons.faces;
-  BOOST_FOREACH(IndexedFace &face, faces) {
-		int i=0;
-		while (i < face.size()) {
+	std::vector<IndexedFace> 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<faces.size();i++) {
-		if (faces[i].size() < 3) {
-			faces.erase(faces.begin() + i);
+	for (size_t i=1;i<cleanfaces.size();i++) {
+		if (cleanfaces[i].size() < 3) {
+			cleanfaces.erase(cleanfaces.begin() + i);
 			i--;
 		}
 	}
 
-	if (faces.size() == 1 && faces[0].size() == 3) {
+	if (cleanfaces.size() == 1 && cleanfaces[0].size() == 3) {
 		// Input polygon has 3 points. shortcut tessellation.
-		triangles.push_back(IndexedTriangle(faces[0][0], faces[0][1], faces[0][2]));
+		//PRINTDB("  tri: %d %d %d", cleanfaces[0][0] % cleanfaces[0][1] % cleanfaces[0][2]);
+		triangles.push_back(IndexedTriangle(cleanfaces[0][0], cleanfaces[0][1], cleanfaces[0][2]));
 		return false;
 	}
 
-	const Vector3f *verts = &polygons.vertices.front();
+	// Build edge dict.
+  // This contains all edges in the original polygon.
+	// To maintain connectivity, all these edges must exist in the output.
+	EdgeDict edges;
+	BOOST_FOREACH(IndexedFace &face, cleanfaces) {
+		edges.add(face);
+	}
+
   TESSreal *normalvec = NULL;
   TESSreal passednormal[3];
   if (normal) {
@@ -97,10 +258,10 @@ bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons,
 	// Since libtess2's indices is based on the running number of points added, we need to map back
 	// to our indices. allindices does the mapping.
 	std::vector<int> 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<inputSize;i++) {
 			if (!vflags[i]) { // vertex missing in output: C)
 				int starti = (i+inputSize-1)%inputSize;
 				int j;
+				PRINTD("   Fanning left-out vertices");
 				for (j = i; j < inputSize && !vflags[j]; j++) {
 					// Create triangle fan from vertex i-1 to the first existing vertex
-					PRINTDB("(%d) (%d) (%d)\n", allindices[starti] % allindices[j] % allindices[((j+1)%inputSize)]);
+					PRINTDB("   (%d) (%d) (%d)\n", allindices[starti] % allindices[j] % allindices[((j+1)%inputSize)]);
 					tri[0] = allindices[starti];
 					tri[1] = allindices[j];
 					tri[2] = allindices[(j+1)%inputSize];
@@ -184,27 +386,7 @@ bool GeometryUtils::tessellatePolygonWithHoles(const IndexedPolygons &polygons,
 				i = j;
 			}
 		}
-	}
-	else {
-		IndexedTriangle tri;
-		for (int t=0;t<numelems;t++) {
-			bool err = false;
-			for (int i=0;i<3;i++) {
-				int vidx = vindices[elements[t*3 + i]];
-				if (vidx == TESS_UNDEF) err = true;
-				else tri[i] = allindices[vidx];
-			}
-			PRINTDB("%d (%d) %d (%d) %d (%d)",
-							elements[t*3 + 0] % allindices[vindices[elements[t*3 + 0]]] %
-							elements[t*3 + 1] % allindices[vindices[elements[t*3 + 1]]] %
-							elements[t*3 + 2] % allindices[vindices[elements[t*3 + 2]]]);
-			// FIXME: We ignore self-intersecting triangles rather than detecting and handling this
-			if (!err) {
-				triangles.push_back(tri);
-			}
-			else PRINT("WARNING: Self-intersecting polygon encountered - ignoring");
-		}
-	}
+#endif
 
   tessDeleteTess(tess);
 
@@ -219,19 +401,18 @@ bool GeometryUtils::tessellatePolygon(const Polygon &polygon, Polygons &triangle
 {
 	bool err = false;
 	Reindexer<Vector3f> uniqueVertices;
-	IndexedPolygons indexedpolygons;
-	indexedpolygons.faces.push_back(IndexedFace());
-	IndexedFace &currface = indexedpolygons.faces.back();
+	std::vector<IndexedFace> indexedfaces;
+	indexedfaces.push_back(IndexedFace());
+	IndexedFace &currface = indexedfaces.back();
 	BOOST_FOREACH (const Vector3d &v, polygon) {
 		int idx = uniqueVertices.lookup(v.cast<float>());
 		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<IndexedTriangle> 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<std::vector<IndexedFace> > &polygons)
+{
+	EdgeDict edges;
+	BOOST_FOREACH(const std::vector<IndexedFace> &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<IndexedTriangle> &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();
+}

src/GeometryUtils.h

@@ -19,9 +19,21 @@ struct IndexedTriangleMesh {
 	std::vector<IndexedTriangle> triangles;
 };
 
+// Indexed polygon mesh, where each polygon can have holes
+struct IndexedPolyMesh {
+	std::vector<Vector3f> vertices;
+	std::vector<std::vector<IndexedFace> > 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<IndexedTriangle> &triangles,
+	bool tessellatePolygonWithHoles(const Vector3f *vertices,
+																	const std::vector<IndexedFace> &faces, 
+																	std::vector<IndexedTriangle> &triangles,
 																	const Vector3f *normal = NULL);
+
+	int findUnconnectedEdges(const std::vector<std::vector<IndexedFace> > &polygons);
+	int findUnconnectedEdges(const std::vector<IndexedTriangle> &triangles);
 }

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<Vector3f> allVertices;
+		std::vector<std::vector<IndexedFace> > polygons;
+
+		CGAL_Nef_polyhedron3::Halffacet_const_iterator hfaceti;
+		CGAL_forall_halffacets(hfaceti, N) {
+			CGAL::Plane_3<CGAL_Kernel3> 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<IndexedFace>());
+			std::vector<IndexedFace> &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<Vector3f>(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<IndexedTriangle> allTriangles;
+		BOOST_FOREACH(const std::vector<IndexedFace> &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<allVertices.size();i++) {
+			std::cerr << verts[i][0] << ", " << verts[i][1] << ", " << verts[i][2] << "\n";
+		}		
+#endif
+
+			/* at this stage, we have a sequence of polygons. the first
+				 is the "outside edge' or 'body' or 'border', and the rest of the
+				 polygons are 'holes' within the first. there are several
+				 options here to get rid of the holes. we choose to go ahead
+				 and let the tessellater deal with the holes, and then
+				 just output the resulting 3d triangles*/
+
+			// We cannot trust the plane from Nef polyhedron to be correct.
+			// Passing an incorrect normal vector can cause a crash in the constrained delaunay triangulator
+      // See http://cgal-discuss.949826.n4.nabble.com/Nef3-Wrong-normal-vector-reported-causes-triangulator-crash-tt4660282.html
+			// CGAL::Vector_3<CGAL_Kernel3> 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<IndexedTriangle> 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<allVertices.size();i++) {
+			std::cerr << verts[i][0] << ", " << verts[i][1] << ", " << verts[i][2] << "\n";
+		}		
+#endif
+
 		return err;
 	}
 #endif