Make sure convex objects remain convex as Nef polyhedrons

src/GeometryEvaluator.cc

@@ -97,7 +97,7 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
 	if (children.size() == 0) return ResultObject();
 
 	if (op == OPENSCAD_HULL) {
-		PolySet *ps = new PolySet(3);
+		PolySet *ps = new PolySet(3, true);
 
 		if (CGALUtils::applyHull(children, *ps)) {
 			return ps;

src/cgalutils.cc

@@ -11,12 +11,20 @@
 #include "cgal.h"
 #include <CGAL/convex_hull_3.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/normal_vector_newell_3.h>
 #include "svg.h"
 #include "Reindexer.h"
 
 #include <map>
 #include <boost/foreach.hpp>
 
+namespace /* anonymous */ {
+	template<typename Result, typename V>
+	Result vector_convert(V const& v) {
+		return Result(CGAL::to_double(v[0]),CGAL::to_double(v[1]),CGAL::to_double(v[2]));
+	}
+}
+
 namespace CGALUtils {
 
 	bool applyHull(const Geometry::ChildList &children, PolySet &result)
@@ -337,6 +345,84 @@ namespace CGALUtils {
 		return result;
 	}
 
+	namespace {
+
+		// lexicographic comparison
+		bool operator < (Vector3d const& a, Vector3d const& b) {
+			for (int i = 0; i < 3; i++) {
+				if (a[i] < b[i]) return true;
+				else if (a[i] == b[i]) continue;
+				return false;
+			}
+			return false;
+		}
+	}
+
+	struct VecPairCompare {
+		bool operator ()(std::pair<Vector3d, Vector3d> const& a,
+						 std::pair<Vector3d, Vector3d> const& b) const {
+			return a.first < b.first || (!(b.first < a.first) && a.second < b.second);
+		}
+	};
+
+
+	bool is_approximately_convex(const PolySet &ps) {
+
+		const double angle_threshold = cos(.5/180*M_PI); // .5°
+
+		typedef CGAL::Simple_cartesian<double> K;
+		typedef K::Vector_3 Vector;
+		typedef K::Point_3 Point;
+		typedef K::Plane_3 Plane;
+
+		// compute edge to face relations and plane equations
+		typedef std::pair<Vector3d,Vector3d> Edge;
+		typedef std::map<Edge, int, VecPairCompare> Edge_to_facet_map;
+		Edge_to_facet_map edge_to_facet_map;
+		std::vector<Plane> facet_planes; facet_planes.reserve(ps.polygons.size());
+
+		for (int i = 0; i < ps.polygons.size(); i++) {
+			size_t N = ps.polygons[i].size();
+			assert(N > 0);
+			std::vector<Point> v(N);
+			for (int j = 0; j < N; j++) {
+				v[j] = vector_convert<Point>(ps.polygons[i][j]);
+				Edge edge(ps.polygons[i][j],ps.polygons[i][(j+1)%N]);
+				if (edge_to_facet_map.count(edge)) return false; // edge already exists: nonmanifold
+				edge_to_facet_map[edge] = i;
+			}
+			Vector normal;
+			CGAL::normal_vector_newell_3(v.begin(), v.end(), normal);
+
+			facet_planes.push_back(Plane(v[0], normal));
+		}
+
+		for (int i = 0; i < ps.polygons.size(); i++) {
+			size_t N = ps.polygons[i].size();
+			for (int j = 0; j < N; j++) {
+				Edge other_edge(ps.polygons[i][(j+1)%N], ps.polygons[i][j]);
+				if (edge_to_facet_map.count(other_edge) == 0) return false;//
+				//Edge_to_facet_map::const_iterator it = edge_to_facet_map.find(other_edge);
+				//if (it == edge_to_facet_map.end()) return false; // not a closed manifold
+				//int other_facet = it->second;
+				int other_facet = edge_to_facet_map[other_edge];
+
+				Point p = vector_convert<Point>(ps.polygons[i][(j+2)%N]);
+
+				if (facet_planes[other_facet].has_on_positive_side(p)) {
+					// Check angle
+					Vector u = facet_planes[other_facet].orthogonal_vector();
+					Vector v = facet_planes[i].orthogonal_vector();
+
+					double cos_angle = u / sqrt(u*u) * v / sqrt(v*v);
+					if (cos_angle < angle_threshold) {
+						return false;
+					}
+				}
+			}
+		}
+		return true;
+	}
 };
 
 template <typename Polyhedron>
@@ -1100,11 +1186,98 @@ void ZRemover::visit( CGAL_Nef_polyhedron3::Halffacet_const_handle hfacet )
 	PRINTD(" <!-- ZRemover Halffacet visit end -->");
 }
 
+namespace /* anonymous */ {
+	// This code is from CGAL/demo/Polyhedron/Scene_nef_polyhedron_item.cpp
+	// quick hacks to convert polyhedra from exact to inexact and vice-versa
+	template <class Polyhedron_input,
+	class Polyhedron_output>
+	struct Copy_polyhedron_to
+	: public CGAL::Modifier_base<typename Polyhedron_output::HalfedgeDS>
+	{
+		Copy_polyhedron_to(const Polyhedron_input& in_poly)
+		: in_poly(in_poly) {}
+
+		void operator()(typename Polyhedron_output::HalfedgeDS& out_hds)
+		{
+			typedef typename Polyhedron_output::HalfedgeDS Output_HDS;
+
+			CGAL::Polyhedron_incremental_builder_3<Output_HDS> builder(out_hds);
+
+			typedef typename Polyhedron_input::Vertex_const_iterator Vertex_const_iterator;
+			typedef typename Polyhedron_input::Facet_const_iterator  Facet_const_iterator;
+			typedef typename Polyhedron_input::Halfedge_around_facet_const_circulator HFCC;
+
+			builder.begin_surface(in_poly.size_of_vertices(),
+								  in_poly.size_of_facets(),
+								  in_poly.size_of_halfedges());
+
+			for(Vertex_const_iterator
+				vi = in_poly.vertices_begin(), end = in_poly.vertices_end();
+				vi != end ; ++vi)
+			{
+				typename Polyhedron_output::Point_3 p(::CGAL::to_double( vi->point().x()),
+													  ::CGAL::to_double( vi->point().y()),
+													  ::CGAL::to_double( vi->point().z()));
+				builder.add_vertex(p);
+			}
+
+			typedef CGAL::Inverse_index<Vertex_const_iterator> Index;
+			Index index( in_poly.vertices_begin(), in_poly.vertices_end());
+
+			for(Facet_const_iterator
+				fi = in_poly.facets_begin(), end = in_poly.facets_end();
+				fi != end; ++fi)
+			{
+				HFCC hc = fi->facet_begin();
+				HFCC hc_end = hc;
+				//     std::size_t n = circulator_size( hc);
+				//     CGAL_assertion( n >= 3);
+				builder.begin_facet ();
+				do {
+					builder.add_vertex_to_facet(index[hc->vertex()]);
+					++hc;
+				} while( hc != hc_end);
+				builder.end_facet();
+			}
+			builder.end_surface();
+		} // end operator()(..)
+	private:
+		const Polyhedron_input& in_poly;
+	}; // end Copy_polyhedron_to<>
+
+	template <class Poly_A, class Poly_B>
+	void copy_to(const Poly_A& poly_a, Poly_B& poly_b)
+	{
+		Copy_polyhedron_to<Poly_A, Poly_B> modifier(poly_a);
+		poly_b.delegate(modifier);
+	}
+}
+
 static CGAL_Nef_polyhedron *createNefPolyhedronFromPolySet(const PolySet &ps)
 {
 	if (ps.isEmpty()) return new CGAL_Nef_polyhedron();
 	assert(ps.getDimension() == 3);
 
+	if (ps.is_convex()) {
+		typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+		// Collect point cloud
+		std::set<K::Point_3> points;
+		for (int i = 0; i < ps.polygons.size(); i++) {
+			for (int j = 0; j < ps.polygons[i].size(); j++) {
+				points.insert(vector_convert<K::Point_3>(ps.polygons[i][j]));
+			}
+		}
+
+		if (points.size() <= 3) return new CGAL_Nef_polyhedron();;
+
+		// Apply hull
+		CGAL::Polyhedron_3<K> r;
+		CGAL::convex_hull_3(points.begin(), points.end(), r);
+		CGAL::Polyhedron_3<CGAL_Kernel3> r_exact;
+		copy_to(r,r_exact);
+		return new CGAL_Nef_polyhedron(new CGAL_Nef_polyhedron3(r_exact));
+	}
+
 	CGAL_Nef_polyhedron3 *N = NULL;
 	bool plane_error = false;
 	CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION);

src/cgalutils.h

@@ -11,6 +11,7 @@ namespace CGALUtils {
 	void applyBinaryOperator(CGAL_Nef_polyhedron &target, const CGAL_Nef_polyhedron &src, OpenSCADOperator op);
 	Polygon2d *project(const CGAL_Nef_polyhedron &N, bool cut);
 	CGAL_Iso_cuboid_3 boundingBox(const CGAL_Nef_polyhedron3 &N);
+	bool is_approximately_convex(const PolySet &ps);
 };
 
 CGAL_Nef_polyhedron *createNefPolyhedronFromGeometry(const class Geometry &geom);

src/polyset.cc

@@ -44,11 +44,11 @@
 
  */
 
-PolySet::PolySet(unsigned int dim) : dim(dim)
+PolySet::PolySet(unsigned int dim, boost::tribool convex) : dim(dim), convex(convex)
 {
 }
 
-PolySet::PolySet(const Polygon2d &origin) : polygon(origin), dim(2)
+PolySet::PolySet(const Polygon2d &origin) : polygon(origin), dim(2), convex(unknown)
 {
 }
 
@@ -140,6 +140,22 @@ void PolySet::transform(const Transform3d &mat)
 	}
 }
 
+namespace CGALUtils {
+	extern bool is_approximately_convex(const PolySet &ps);
+}
+
+bool PolySet::is_convex() const {
+	if (convex)  return true;
+	if (!convex) return false;
+
+#ifdef ENABLE_CGAL
+	convex = CGALUtils::is_approximately_convex(*this);
+	return convex;
+#else
+	return false;
+#endif
+}
+
 void PolySet::resize(Vector3d newsize, const Eigen::Matrix<bool,3,1> &autosize)
 {
 	BoundingBox bbox = this->getBoundingBox();

src/polyset.h

@@ -8,13 +8,16 @@
 #include <vector>
 #include <string>
 
+#include <boost/logic/tribool.hpp>
+BOOST_TRIBOOL_THIRD_STATE(unknown)
+
 class PolySet : public Geometry
 {
 public:
 	typedef std::vector<Vector3d> Polygon;
 	std::vector<Polygon> polygons;
 
-	PolySet(unsigned int dim);
+	PolySet(unsigned int dim, boost::tribool convex = unknown);
 	PolySet(const Polygon2d &origin);
 	virtual ~PolySet();
 
@@ -38,7 +41,10 @@ public:
 	void transform(const Transform3d &mat);
 	void resize(Vector3d newsize, const Eigen::Matrix<bool,3,1> &autosize);
 
+	bool is_convex() const;
+
 private:
-  Polygon2d polygon;
+	Polygon2d polygon;
 	unsigned int dim;
+	mutable boost::tribool convex;
 };

src/primitives.cc

@@ -305,7 +305,7 @@ Geometry *PrimitiveNode::createGeometry() const
 
 	switch (this->type) {
 	case CUBE: {
-		PolySet *p = new PolySet(3);
+		PolySet *p = new PolySet(3,true);
 		g = p;
 		if (this->x > 0 && this->y > 0 && this->z > 0 &&
 			!isinf(this->x) > 0 && !isinf(this->y) > 0 && !isinf(this->z) > 0) {
@@ -363,7 +363,7 @@ Geometry *PrimitiveNode::createGeometry() const
 	}
 		break;
 	case SPHERE: {
-		PolySet *p = new PolySet(3);
+		PolySet *p = new PolySet(3,true);
 		g = p;
 		if (this->r1 > 0 && !isinf(this->r1)) {
 			struct ring_s {
@@ -438,7 +438,7 @@ Geometry *PrimitiveNode::createGeometry() const
 	}
 		break;
 	case CYLINDER: {
-		PolySet *p = new PolySet(3);
+		PolySet *p = new PolySet(3,true);
 		g = p;
 		if (this->h > 0 && !isinf(this->h) &&
 				this->r1 >=0 && this->r2 >= 0 && (this->r1 > 0 || this->r2 > 0) &&