diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/CGALEvaluator.cc | 26 | ||||
-rw-r--r-- | src/CGAL_Nef_polyhedron.cc | 5 | ||||
-rw-r--r-- | src/cgalutils.cc | 27 | ||||
-rw-r--r-- | src/cgalutils.h | 5 | ||||
-rw-r--r-- | src/dxftess-cgal.cc | 175 | ||||
-rw-r--r-- | src/dxftess.h | 1 | ||||
-rw-r--r-- | src/import.cc | 4 | ||||
-rw-r--r-- | src/polyset.h | 2 | ||||
-rw-r--r-- | src/primitives.cc | 18 |
9 files changed, 230 insertions, 33 deletions
diff --git a/src/CGALEvaluator.cc b/src/CGALEvaluator.cc index 86118d7..20c5d5e 100644 --- a/src/CGALEvaluator.cc +++ b/src/CGALEvaluator.cc @@ -680,16 +680,32 @@ CGAL_Nef_polyhedron CGALEvaluator::evaluateCGALMesh(const PolySet &ps) else // not (this->is2d) { CGAL_Nef_polyhedron3 *N = NULL; + bool plane_error = false; CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); try { - // FIXME: Are we leaking memory for the CGAL_Polyhedron object? - CGAL_Polyhedron *P = createPolyhedronFromPolySet(ps); - if (P) { - N = new CGAL_Nef_polyhedron3(*P); + CGAL_Polyhedron P; + bool err = createPolyhedronFromPolySet(ps,P); + if (!err) N = new CGAL_Nef_polyhedron3(P); + } + catch (const CGAL::Assertion_exception &e) { + if (std::string(e.what()).find("Plane_constructor")!=std::string::npos) { + if (std::string(e.what()).find("has_on")!=std::string::npos) { + PRINT("PolySet has nonplanar faces. Attempting alternate construction"); + plane_error=true; + } + } else { + PRINTB("CGAL error in CGAL_Nef_polyhedron3(): %s", e.what()); } } + if (plane_error) try { + PolySet ps2; + CGAL_Polyhedron P; + tessellate_3d_faces( ps, ps2 ); + bool err = createPolyhedronFromPolySet(ps2,P); + if (!err) N = new CGAL_Nef_polyhedron3(P); + } catch (const CGAL::Assertion_exception &e) { - PRINTB("CGAL error in CGAL_Nef_polyhedron3(): %s", e.what()); + PRINTB("Alternate construction failed. CGAL error in CGAL_Nef_polyhedron3(): %s", e.what()); } CGAL::set_error_behaviour(old_behaviour); return CGAL_Nef_polyhedron(N); diff --git a/src/CGAL_Nef_polyhedron.cc b/src/CGAL_Nef_polyhedron.cc index 49b9a53..ea9accd 100644 --- a/src/CGAL_Nef_polyhedron.cc +++ b/src/CGAL_Nef_polyhedron.cc @@ -96,6 +96,7 @@ PolySet *CGAL_Nef_polyhedron::convertToPolyset() } else if (this->dim == 3) { CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); + ps = new PolySet(); bool err = true; std::string errmsg(""); CGAL_Polyhedron P; @@ -107,11 +108,11 @@ PolySet *CGAL_Nef_polyhedron::convertToPolyset() err = true; errmsg = std::string(e.what()); } + if (!err) err = createPolySetFromPolyhedron(P, *ps); if (err) { PRINT("ERROR: CGAL NefPolyhedron->Polyhedron conversion failed."); if (errmsg!="") PRINTB("ERROR: %s",errmsg); - } else { - ps = createPolySetFromPolyhedron(P); + delete ps; ps = NULL; } CGAL::set_error_behaviour(old_behaviour); } 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 d25fd4c..8f7e4dd 100644 --- a/src/cgalutils.h +++ b/src/cgalutils.h @@ -2,8 +2,9 @@ #define CGALUTILS_H_ #include <cgal.h> -class PolySet *createPolySetFromPolyhedron(const CGAL_Polyhedron &p); -CGAL_Polyhedron *createPolyhedronFromPolySet(const class PolySet &ps); +#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 16eaf9f..14f1204 100644 --- a/src/dxftess-cgal.cc +++ b/src/dxftess-cgal.cc @@ -335,3 +335,178 @@ void dxf_tesselate(PolySet *ps, DxfData &dxf, double rot, Vector2d scale, bool u dxf.paths[path[2]].is_inner = !up; } } + + +/* Tessellation of 3d PolySet faces + +This code is for tessellating the faces of a 3d PolySet, assuming that +the faces are near-planar polygons. + +We do the tessellation by projecting each polygon of the Polyset onto a +2-d plane, then running a 2d tessellation algorithm on the projected 2d +polygon. Then we project each of the newly generated 2d 'tiles' (the +polygons used for tessellation, typically triangles) back up into 3d +space. + +(in reality as of writing, we dont need to do a back-projection from 2d->3d +because the algorithm we are using doesn't create any new points, and we can +just use a 'map' to associate 3d points with 2d points). + +The code assumes the input polygons are simple, non-intersecting, without +holes, without duplicate input points, and with proper orientation. + +The purpose of this code is originally to fix github issue 349. Our CGAL +kernel does not accept polygons for Nef_Polyhedron_3 if each of the +points is not exactly coplanar. "Near-planar" or "Almost planar" polygons +often occur due to rounding issues on, for example, polyhedron() input. +By tessellating the 3d polygon into individual smaller tiles that +are perfectly coplanar (triangles, for example), we can get CGAL to accept +the polyhedron() input. +*/ + +typedef enum { XYPLANE, YZPLANE, XZPLANE, NONE } projection_t; + +// this is how we make 3d points appear as though they were 2d points to +//the tessellation algorithm. +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() ); } + +/* Find a 'good' 2d projection for a given 3d polygon. the XY, YZ, or XZ +plane. This is needed because near-planar polygons in 3d can have 'bad' +projections into 2d. For example if the square 0,0,0 0,1,0 0,1,1 0,0,1 +is projected onto the XY plane you will not get a polygon, you wil get +a skinny line thing. It's better to project that square onto the yz +plane.*/ +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 XY, YZ, or XZ + // 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; +} + +/* triangulate the given 3d polygon using CGAL's 2d Constrained Delaunay +algorithm. Project the polygon's points into 2d using the given projection +before performing the triangulation. This code assumes input polygon is +simple, no holes, no self-intersections, no duplicate points, and is +properly oriented. output is a sequence of 3d triangles. */ +bool triangulate_polygon( const PolySet::Polygon &pgon, std::vector<PolySet::Polygon> &triangles, projection_t projection ) +{ + bool err = false; + CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); + try { + CDT cdt; + std::vector<Vertex_handle> vhandles; + std::map<CDTPoint,Vector3d> vertmap; + CGAL::Orientation original_orientation; + std::vector<CDTPoint> orienpgon; + 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); + orienpgon.push_back( cdtpoint ); + } + original_orientation = CGAL::orientation_2( orienpgon.begin(),orienpgon.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 ); + } + } + } catch (const CGAL::Assertion_exception &e) { + PRINTB("CGAL error in dxftess triangulate_polygon: %s", e.what()); + err = true; + } + CGAL::set_error_behaviour(old_behaviour); + return err; +} + +/* Given a 3d PolySet with 'near planar' polygonal faces, Tessellate the +faces. As of writing, our only tessellation method is Triangulation +using CGAL's Constrained Delaunay algorithm. This code assumes the input +polyset has simple polygon faces with no holes, no self intersections, no +duplicate points, and proper orientation. */ +void tessellate_3d_faces( const PolySet &inps, PolySet &outps ) { + for (size_t i = 0; i < inps.polygons.size(); i++) { + const PolySet::Polygon pgon = inps.polygons[i]; + if (pgon.size()<3) { + PRINT("WARNING: PolySet has polygon with <3 points"); + continue; + } + projection_t goodproj = find_good_projection( pgon ); + if (goodproj==NONE) { + PRINT("WARNING: PolySet has degenerate polygon"); + continue; + } + std::vector<PolySet::Polygon> triangles; + bool err = triangulate_polygon( pgon, triangles, goodproj ); + if (!err) 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()); + } + } +} + +// End of PolySet face tessellation code + diff --git a/src/dxftess.h b/src/dxftess.h index f0f27b5..d596f0b 100644 --- a/src/dxftess.h +++ b/src/dxftess.h @@ -7,5 +7,6 @@ class DxfData; class PolySet; void dxf_tesselate(PolySet *ps, DxfData &dxf, double rot, Vector2d scale, bool up, bool do_triangle_splitting, double h); void dxf_border_to_ps(PolySet *ps, const DxfData &dxf); +void tessellate_3d_faces( const PolySet &inps, PolySet &outps ); #endif diff --git a/src/import.cc b/src/import.cc index 3897331..3c2cce1 100644 --- a/src/import.cc +++ b/src/import.cc @@ -284,7 +284,9 @@ PolySet *ImportNode::evaluate_polyset(class PolySetEvaluator *) const file >> poly; file.close(); - p = createPolySetFromPolyhedron(poly); + p = new PolySet(); + bool err = createPolySetFromPolyhedron(poly, *p); + if (err) delete p; } #else PRINT("WARNING: OFF import requires CGAL."); diff --git a/src/polyset.h b/src/polyset.h index 6626f79..2ff817a 100644 --- a/src/polyset.h +++ b/src/polyset.h @@ -26,7 +26,7 @@ public: void append_vertex(double x, double y, double z = 0.0); void insert_vertex(double x, double y, double z = 0.0); size_t memsize() const; - + BoundingBox getBoundingBox() const; #define CSGMODE_DIFFERENCE_FLAG 0x10 diff --git a/src/primitives.cc b/src/primitives.cc index a587d43..a76637d 100644 --- a/src/primitives.cc +++ b/src/primitives.cc @@ -103,7 +103,7 @@ public: double fn, fs, fa; primitive_type_e type; int convexity; - Value points, paths, triangles; + Value points, paths, faces; virtual PolySet *evaluate_polyset(class PolySetEvaluator *) const; }; @@ -156,7 +156,7 @@ AbstractNode *PrimitiveModule::instantiate(const Context *ctx, const ModuleInsta args += Assignment("h", NULL), Assignment("r1", NULL), Assignment("r2", NULL), Assignment("center", NULL); break; case POLYHEDRON: - args += Assignment("points", NULL), Assignment("triangles", NULL), Assignment("convexity", NULL); + args += Assignment("points", NULL), Assignment("triangles", NULL), Assignment("convexity", NULL), Assignment("faces", NULL); break; case SQUARE: args += Assignment("size", NULL), Assignment("center", NULL); @@ -235,7 +235,11 @@ AbstractNode *PrimitiveModule::instantiate(const Context *ctx, const ModuleInsta if (type == POLYHEDRON) { node->points = c.lookup_variable("points"); - node->triangles = c.lookup_variable("triangles"); + node->faces = c.lookup_variable("faces"); + if (node->faces.type() == Value::UNDEFINED) { + // backwards compatable + node->faces = c.lookup_variable("triangles"); + } } if (type == SQUARE) { @@ -480,11 +484,11 @@ sphere_next_r2: if (this->type == POLYHEDRON) { p->convexity = this->convexity; - for (size_t i=0; i<this->triangles.toVector().size(); i++) + for (size_t i=0; i<this->faces.toVector().size(); i++) { p->append_poly(); - for (size_t j=0; j<this->triangles.toVector()[i].toVector().size(); j++) { - size_t pt = this->triangles.toVector()[i].toVector()[j].toDouble(); + for (size_t j=0; j<this->faces.toVector()[i].toVector().size(); j++) { + size_t pt = this->faces.toVector()[i].toVector()[j].toDouble(); if (pt < this->points.toVector().size()) { double px, py, pz; if (this->points.toVector()[pt].getVec3(px, py, pz)) @@ -610,7 +614,7 @@ std::string PrimitiveNode::toString() const break; case POLYHEDRON: stream << "(points = " << this->points - << ", triangles = " << this->triangles + << ", faces = " << this->faces << ", convexity = " << this->convexity << ")"; break; case SQUARE: |