From 8c977a45e579a7b72d44347ea2916105bf6b9cd4 Mon Sep 17 00:00:00 2001 From: Oskar Linde Date: Sun, 1 Jun 2014 00:35:05 +0200 Subject: [PATCH] Make sure convex objects remain convex as Nef polyhedrons --- src/GeometryEvaluator.cc | 2 +- src/cgalutils.cc | 173 +++++++++++++++++++++++++++++++++++++++ src/cgalutils.h | 1 + src/polyset.cc | 20 ++++- src/polyset.h | 10 ++- src/primitives.cc | 6 +- 6 files changed, 204 insertions(+), 8 deletions(-) diff --git a/src/GeometryEvaluator.cc b/src/GeometryEvaluator.cc index 4af9d0c0..df6ff0fe 100644 --- a/src/GeometryEvaluator.cc +++ b/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; diff --git a/src/cgalutils.cc b/src/cgalutils.cc index 13b6a66e..248992d5 100644 --- a/src/cgalutils.cc +++ b/src/cgalutils.cc @@ -11,12 +11,20 @@ #include "cgal.h" #include #include +#include #include "svg.h" #include "Reindexer.h" #include #include +namespace /* anonymous */ { + template + 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 const& a, + std::pair 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 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 Edge; + typedef std::map Edge_to_facet_map; + Edge_to_facet_map edge_to_facet_map; + std::vector 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 v(N); + for (int j = 0; j < N; j++) { + v[j] = vector_convert(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(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 @@ -1100,11 +1186,98 @@ void ZRemover::visit( CGAL_Nef_polyhedron3::Halffacet_const_handle hfacet ) PRINTD(" "); } +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 + struct Copy_polyhedron_to + : public CGAL::Modifier_base + { + 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 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 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 + void copy_to(const Poly_A& poly_a, Poly_B& poly_b) + { + Copy_polyhedron_to 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 points; + for (int i = 0; i < ps.polygons.size(); i++) { + for (int j = 0; j < ps.polygons[i].size(); j++) { + points.insert(vector_convert(ps.polygons[i][j])); + } + } + + if (points.size() <= 3) return new CGAL_Nef_polyhedron();; + + // Apply hull + CGAL::Polyhedron_3 r; + CGAL::convex_hull_3(points.begin(), points.end(), r); + CGAL::Polyhedron_3 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); diff --git a/src/cgalutils.h b/src/cgalutils.h index 29295cad..432b1ce5 100644 --- a/src/cgalutils.h +++ b/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); diff --git a/src/polyset.cc b/src/polyset.cc index 689232c4..0d811e4a 100644 --- a/src/polyset.cc +++ b/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 &autosize) { BoundingBox bbox = this->getBoundingBox(); diff --git a/src/polyset.h b/src/polyset.h index c481d5aa..852084ef 100644 --- a/src/polyset.h +++ b/src/polyset.h @@ -8,13 +8,16 @@ #include #include +#include +BOOST_TRIBOOL_THIRD_STATE(unknown) + class PolySet : public Geometry { public: typedef std::vector Polygon; std::vector 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 &autosize); + bool is_convex() const; + private: - Polygon2d polygon; + Polygon2d polygon; unsigned int dim; + mutable boost::tribool convex; }; diff --git a/src/primitives.cc b/src/primitives.cc index 401e3304..b0ad1be6 100644 --- a/src/primitives.cc +++ b/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) &&