From a49c32bee056f318b0910934f7b8e433c7aef189 Mon Sep 17 00:00:00 2001 From: Marius Kintel Date: Sun, 15 Dec 2013 21:53:24 -0500 Subject: [PATCH] Fixes remaining issues after merging #574 --- src/cgalutils.cc | 30 ++++-- src/dxftess-cgal.cc | 175 ----------------------------------- src/dxftess.h | 1 - src/import.cc | 3 +- src/polyset-utils.cc | 213 +++++++++++++++++++++++++++++++++++++++++++ src/polyset-utils.h | 1 + tests/CMakeLists.txt | 2 +- 7 files changed, 240 insertions(+), 185 deletions(-) diff --git a/src/cgalutils.cc b/src/cgalutils.cc index 40551acb..d1fa323f 100644 --- a/src/cgalutils.cc +++ b/src/cgalutils.cc @@ -404,7 +404,7 @@ static CGAL_Nef_polyhedron *createNefPolyhedronFromPolySet(const PolySet &ps) std::list< point_list_t > pdata_point_lists; std::list < std::pair < point_list_it, point_list_it > > pdata; Grid2d grid(GRID_COARSE); - + for (int i = 0; i < ps.polygons.size(); i++) { pdata_point_lists.push_back(point_list_t()); for (int j = 0; j < ps.polygons[i].size(); j++) { @@ -643,16 +643,32 @@ static CGAL_Nef_polyhedron *createNefPolyhedronFromPolySet(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) { - PRINTB("CGAL error in CGALUtils::createNefPolyhedronFromPolySet(): %s", e.what()); + 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; + PolysetUtils::tessellate_faces(ps, ps2); + bool err = createPolyhedronFromPolySet(ps2,P); + if (!err) N = new CGAL_Nef_polyhedron3(P); + } + catch (const CGAL::Assertion_exception &e) { + PRINTB("Alternate construction failed. CGAL error in CGAL_Nef_polyhedron3(): %s", e.what()); } CGAL::set_error_behaviour(old_behaviour); return new CGAL_Nef_polyhedron(N); diff --git a/src/dxftess-cgal.cc b/src/dxftess-cgal.cc index 5f95e6a0..16eaf9fb 100644 --- a/src/dxftess-cgal.cc +++ b/src/dxftess-cgal.cc @@ -335,178 +335,3 @@ 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 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 &triangles, projection_t projection ) -{ - bool err = false; - CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); - try { - CDT cdt; - std::vector vhandles; - std::map vertmap; - CGAL::Orientation original_orientation; - std::vector 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 list_of_seeds; - CGAL::refine_Delaunay_mesh_2_without_edge_refinement(cdt, - list_of_seeds.begin(), list_of_seeds.end(), DummyCriteria()); - - 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::Failure_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 triangles; - bool err = triangulate_polygon( pgon, triangles, goodproj ); - if (!err) for (size_t j=0;j> poly; file.close(); - p = new PolySet(); + PolySet *p = new PolySet(); bool err = createPolySetFromPolyhedron(poly, *p); if (err) delete p; + else g = p; } #else PRINT("WARNING: OFF import requires CGAL."); diff --git a/src/polyset-utils.cc b/src/polyset-utils.cc index 78b56e07..7a01bf04 100644 --- a/src/polyset-utils.cc +++ b/src/polyset-utils.cc @@ -1,6 +1,49 @@ #include "polyset-utils.h" #include "polyset.h" #include "Polygon2d.h" +#include "printutils.h" +#include "cgal.h" + +#ifdef NDEBUG +#define PREV_NDEBUG NDEBUG +#undef NDEBUG +#endif +#include +#include +#include +#include +#include +#include +#include +#ifdef PREV_NDEBUG +#define NDEBUG PREV_NDEBUG +#endif + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Triangulation_vertex_base_2 Vb; +typedef CGAL::Delaunay_mesh_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +//typedef CGAL::Delaunay_mesh_criteria_2 Criteria; + +typedef CDT::Vertex_handle Vertex_handle; +typedef CDT::Point CDTPoint; + +template class DummyCriteria { +public: + typedef double Quality; + class Is_bad { + public: + CGAL::Mesh_2::Face_badness operator()(const Quality) const { + return CGAL::Mesh_2::NOT_BAD; + } + CGAL::Mesh_2::Face_badness operator()(const typename T::Face_handle&, Quality&q) const { + q = 1; + return CGAL::Mesh_2::NOT_BAD; + } + }; + Is_bad is_bad_object() const { return Is_bad(); } +}; #include @@ -22,5 +65,175 @@ namespace PolysetUtils { return poly; } +/* 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 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 &triangles, projection_t projection ) + { + bool err = false; + CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); + try { + CDT cdt; + std::vector vhandles; + std::map vertmap; + CGAL::Orientation original_orientation; + std::vector 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 list_of_seeds; + CGAL::refine_Delaunay_mesh_2_without_edge_refinement(cdt, + list_of_seeds.begin(), list_of_seeds.end(), DummyCriteria()); + + 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_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 triangles; + bool err = triangulate_polygon( pgon, triangles, goodproj ); + if (!err) for (size_t j=0;j