triangulation of near-planar PolySet faces

4 files changed, 127 insertions, 26 deletions

diff --git a/src/CGAL_Nef_polyhedron.cc b/src/CGAL_Nef_polyhedron.cc
index 2d5bba5..8b54eba 100644
--- a/src/CGAL_Nef_polyhedron.cc
+++ b/src/CGAL_Nef_polyhedron.cc
@@ -100,7 +100,7 @@ PolySet *CGAL_Nef_polyhedron::convertToPolyset()
 			ps = new PolySet();
 			CGAL_Polyhedron P;
 			this->p3->convert_to_Polyhedron(P);
-			bool err = createPolySetFromPolyhedron(P, ps);
+			bool err = createPolySetFromPolyhedron(P, *ps);
 			if (err) delete ps;
 		}
 		catch (const CGAL::Precondition_exception &e) {
diff --git a/src/cgalutils.cc b/src/cgalutils.cc
index bb46f1c..12e743d 100644
--- a/src/cgalutils.cc
+++ b/src/cgalutils.cc
@@ -8,10 +8,9 @@
 
 #include <map>
 
-PolySet *createPolySetFromPolyhedron(const CGAL_Polyhedron &p)
+bool createPolySetFromPolyhedron(const CGAL_Polyhedron &p, PolySet &ps)
 {
-	PolySet *ps = new PolySet();
-		
+	bool err = false;
 	typedef CGAL_Polyhedron::Vertex                                 Vertex;
 	typedef CGAL_Polyhedron::Vertex_const_iterator                  VCI;
 	typedef CGAL_Polyhedron::Facet_const_iterator                   FCI;
@@ -35,13 +34,13 @@ PolySet *createPolySetFromPolyhedron(const CGAL_Polyhedron &p)
 			double x3 = CGAL::to_double(v3.point().x());
 			double y3 = CGAL::to_double(v3.point().y());
 			double z3 = CGAL::to_double(v3.point().z());
-			ps->append_poly();
-			ps->append_vertex(x1, y1, z1);
-			ps->append_vertex(x2, y2, z2);
-			ps->append_vertex(x3, y3, z3);
+			ps.append_poly();
+			ps.append_vertex(x1, y1, z1);
+			ps.append_vertex(x2, y2, z2);
+			ps.append_vertex(x3, y3, z3);
 		} while (hc != hc_end);
 	}
-	return ps;
+	return err;
 }
 
 #undef GEN_SURFACE_DEBUG
@@ -127,22 +126,20 @@ public:
 	}
 };
 
-CGAL_Polyhedron *createPolyhedronFromPolySet(const PolySet &ps)
+bool createPolyhedronFromPolySet(const PolySet &ps, CGAL_Polyhedron &p)
 {
-	CGAL_Polyhedron *P = NULL;
+	bool err = false;
 	CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION);
 	try {
-		P = new CGAL_Polyhedron;
 		CGAL_Build_PolySet builder(ps);
-		P->delegate(builder);
+		p.delegate(builder);
 	}
 	catch (const CGAL::Assertion_exception &e) {
 		PRINTB("CGAL error in CGAL_Build_PolySet: %s", e.what());
-		delete P;
-		P = NULL;
+		err = true;
 	}
 	CGAL::set_error_behaviour(old_behaviour);
-	return P;
+	return err;
 }
 
 CGAL_Iso_cuboid_3 bounding_box( const CGAL_Nef_polyhedron3 &N )
diff --git a/src/cgalutils.h b/src/cgalutils.h
index 0e7fbb5..8f7e4dd 100644
--- a/src/cgalutils.h
+++ b/src/cgalutils.h
@@ -2,8 +2,9 @@
 #define CGALUTILS_H_
 
 #include <cgal.h>
-bool createPolySetFromPolyhedron(const CGAL_Polyhedron &p, class PolySet &ps);
-bool createPolyhedronFromPolySet(const class PolySet &ps, CGAL_Polyhedron &p);
+#include "polyset.h"
+bool createPolySetFromPolyhedron(const CGAL_Polyhedron &p, PolySet &ps);
+bool createPolyhedronFromPolySet(const PolySet &ps, CGAL_Polyhedron &p);
 CGAL_Iso_cuboid_3 bounding_box( const CGAL_Nef_polyhedron3 &N );
 CGAL_Iso_rectangle_2e bounding_box( const CGAL_Nef_polyhedron2 &N );
 
diff --git a/src/dxftess-cgal.cc b/src/dxftess-cgal.cc
index d6f592f..6899598 100644
--- a/src/dxftess-cgal.cc
+++ b/src/dxftess-cgal.cc
@@ -336,19 +336,122 @@ void dxf_tesselate(PolySet *ps, DxfData &dxf, double rot, Vector2d scale, bool u
 	}
 }
 
-void triangulate_polygon( const PolySet::Polygon &pgon, std::vector<PolySet::Polygon> &triangles )
-{
+typedef enum { XYPLANE, YZPLANE, XZPLANE, NONE } projection_t;
+
+Vector2d get_projected_point( Vector3d v, projection_t projection ) {
+	Vector2d v2(0,0);
+	if (projection==XYPLANE) { v2.x() = v.x(); v2.y() = v.y(); }
+	else if (projection==XZPLANE) { v2.x() = v.x(); v2.y() = v.z(); }
+	else if (projection==YZPLANE) { v2.x() = v.y(); v2.y() = v.z(); }
+	return v2;
+}
 
+CGAL_Point_3 cgp( Vector3d v ) { return CGAL_Point_3( v.x(), v.y(), v.z() ); }
+/* near-planar polygons in 3d can be projected into 2d, but you have to
+be careful. if you project, say, 0,0,0 0,1,0 0,1,1 0,0,1 onto the XY plane
+you will not get a polygon, you will get a skinny line thing. */
+projection_t find_good_projection( PolySet::Polygon pgon ) {
+	// step 1 - find 3 non-collinear points in the input
+	if (pgon.size()<3) return NONE;
+	Vector3d v1,v2,v3;
+	v1 = v2 = v3 = pgon[0];
+	for (size_t i=0;i<pgon.size();i++) {
+		if (pgon[i]!=v1) { v2=pgon[i]; break; }
+	}
+	if (v1==v2) return NONE;
+	for (size_t i=0;i<pgon.size();i++) {
+		if (!CGAL::collinear( cgp(v1), cgp(v2), cgp(pgon[i]) )) {
+			v3=pgon[i]; break;
+		}
+	}
+	if (CGAL::collinear( cgp(v1), cgp(v2), cgp(v3) ) ) return NONE;
+	// step 2 - find which direction is best for projection. planes use
+	// the equation ax+by+cz+d = 0. a,b, and c determine the direction the
+	// plane is in. we want to find which projection of the 'normal vector'
+	// would make the smallest shadow if projected onto the given plane.
+	// 'quadrance' (distance squared) can tell this w/o using sqrt.
+	CGAL::Plane_3<CGAL_Kernel3> pl( cgp(v1), cgp(v2), cgp(v3) );
+	NT3 qxy = pl.a()*pl.a()+pl.b()*pl.b();
+        NT3 qyz = pl.b()*pl.b()+pl.c()*pl.c();
+        NT3 qxz = pl.c()*pl.c()+pl.a()*pl.a();
+	NT3 min = std::min(qxy,std::min(qyz,qxz));
+	if (min==qxy) return XYPLANE;
+        else if (min==qyz) return YZPLANE;
+        return XZPLANE;
 }
 
-/* given a 3d PolySet with 'near planar' faces, triangulate the faces
-so CGAL Nef Polyhedron will accept them. */
+void triangulate_polygon( const PolySet::Polygon &pgon, std::vector<PolySet::Polygon> &triangles, projection_t projection )
+{
+	CDT cdt;
+	std::vector<Vertex_handle> vhandles;
+	std::map<CDTPoint,Vector3d> vertmap;
+	CGAL::Orientation original_orientation;
+	std::vector<CDTPoint> orpgon;
+	for (size_t i = 0; i < pgon.size(); i++) {
+		Vector3d v3 = pgon.at(i);
+		Vector2d v2 = get_projected_point( v3, projection );
+		CDTPoint cdtpoint = CDTPoint(v2.x(),v2.y());
+		vertmap[ cdtpoint ] = v3;
+		Vertex_handle vh = cdt.insert( cdtpoint );
+		vhandles.push_back(vh);
+		orpgon.push_back( cdtpoint );
+	}
+	original_orientation = CGAL::orientation_2( orpgon.begin(),orpgon.end() );
+	for (size_t i = 0; i < vhandles.size(); i++ ) {
+		int vindex1 = (i+0);
+		int vindex2 = (i+1)%vhandles.size();
+		cdt.insert_constraint( vhandles[vindex1], vhandles[vindex2] );
+	}
+	std::list<CDTPoint> list_of_seeds;
+	CGAL::refine_Delaunay_mesh_2_without_edge_refinement(cdt,
+	  list_of_seeds.begin(), list_of_seeds.end(), DummyCriteria<CDT>());
+
+	CDT::Finite_faces_iterator fit;
+        for( fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); fit++ )
+        {
+                if(fit->is_in_domain()) {
+			CDTPoint p1 = cdt.triangle( fit )[0];
+			CDTPoint p2 = cdt.triangle( fit )[1];
+			CDTPoint p3 = cdt.triangle( fit )[2];
+			Vector3d v1 = vertmap[p1];
+			Vector3d v2 = vertmap[p2];
+			Vector3d v3 = vertmap[p3];
+			PolySet::Polygon pgon;
+			if (CGAL::orientation(p1,p2,p3)==original_orientation) {
+				pgon.push_back(v1);
+				pgon.push_back(v2);
+				pgon.push_back(v3);
+			} else {
+				pgon.push_back(v3);
+				pgon.push_back(v2);
+				pgon.push_back(v1);
+			}
+			triangles.push_back( pgon );
+                }
+        }
+}
+
+/* given a 3d PolySet with 'near planar' faces, triangulate the faces.
+See Issue 340... this is so CGAL Nef Polyhedron will accept them.
+This code assumes the input polyset has simple polygon faces with no holes. */
 void tessellate_3d_faces( const PolySet &inps, PolySet &outps ) {
+	PRINTB("tess 3d %i",inps.polygons.size());
+	PRINTB("%s < input ps",inps.dump());
         for (size_t i = 0; i < inps.polygons.size(); i++) {
-                const PolySet::Polygon *pgon = &inps.polygons[i];
-                for (size_t j = 0; j < pgon->size(); j++) {
-                        Vector3d v = pgon->at(j);
-                }
+                const PolySet::Polygon pgon = inps.polygons[i];
+		if (pgon.size()<3) continue;
+		std::vector<PolySet::Polygon> triangles;
+		projection_t projection = find_good_projection( pgon );
+		triangulate_polygon( pgon, triangles, projection );
+		for (size_t j=0;j<triangles.size();j++) {
+			PolySet::Polygon t = triangles[j];
+			outps.append_poly();
+			outps.append_vertex(t[0].x(),t[0].y(),t[0].z());
+			outps.append_vertex(t[1].x(),t[1].y(),t[1].z());
+			outps.append_vertex(t[2].x(),t[2].y(),t[2].z());
+		}
         }
+	PRINTB("tess 3d done %i",outps.polygons.size());
+	PRINTB("%s < output ps",outps.dump());
 }