diff --git a/src/cgalutils.cc b/src/cgalutils.cc index 5f28e14f..92a779e3 100644 --- a/src/cgalutils.cc +++ b/src/cgalutils.cc @@ -260,6 +260,494 @@ bool createPolySetFromPolyhedron(const CGAL_Polyhedron &p, PolySet &ps) return err; } +/////// Tessellation begin + +typedef CGAL::Plane_3 CGAL_Plane_3; + + +/* + +This is our custom tessellator of Nef Polyhedron faces. The problem with +Nef faces is that sometimes the 'default' tessellator of Nef Polyhedron +doesnt work. This is particularly true with situations where the polygon +face is not, actually, 'simple', according to CGAL itself. This can +occur on a bad quality STL import but also for other reasons. The +resulting Nef face will appear to the average human eye as an ordinary, +simple polygon... but in reality it has multiple edges that are +slightly-out-of-alignment and sometimes they backtrack on themselves. + +When the triangulator is fed a polygon with self-intersecting edges, +it's default behavior is to throw an exception. The other terminology +for this is to say that the 'constraints' in the triangulation are +'intersecting'. The 'constraints' represent the edges of the polygon. +The 'triangulation' is the covering of all the polygon points with +triangles. + +How do we allow interseting constraints during triangulation? We use an +'Itag' for the triangulation, per the CGAL docs. This allows the +triangulator to run without throwing an exception when it encounters +self-intersecting polygon edges. The trick here is that when it finds +an intersection, it actually creates a new point. + +The triangulator creates new points in 2d, but they aren't matched to +any 3d points on our 3d polygon plane. (The plane of the Nef face). How +to fix this problem? We actually 'project back up' or 'lift' into the 3d +plane from the 2d point. This is handled in the 'deproject()' function. + +There is also the issue of the Simplicity of Nef Polyhedron face +polygons. They are often not simple. The intersecting-constraints +Triangulation can triangulate non-simple polygons, but of course it's +result is also non-simple. This means that CGAL functions like +orientation_2() and bounded_side() simply will not work on the resulting +polygons because they all require input polygons to pass the +'is_simple2()' test. We have to use alternatives in order to create our +triangles. + +There is also the question of which underlying number type to use. Some +of the CGAL functions simply dont guarantee good results with a type +like double. Although much the math here is somewhat simple, like +line-line intersection, and involves only simple algebra, the +approximations required when using floating-point types can cause the +answers to be wrong. For example questions like 'is a point inside a +triangle' do not have good answers under floating-point systems where a +line may have a slope that is not expressible exactly as a floating +point number. There are ways to deal with floating point inaccuracy but +it is much, much simpler to use Rational numbers, although potentially +much slower in many cases. + +*/ + +#include +#include + +typedef CGAL_Kernel3 Kernel; +typedef typename CGAL::Triangulation_vertex_base_2 Vb; +//typedef typename CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Delaunay_mesh_face_base_2 Fb; +typedef typename CGAL::Triangulation_data_structure_2 TDS; +typedef CGAL::Exact_intersections_tag ITAG; +typedef typename CGAL::Constrained_Delaunay_triangulation_2 CDT; +//typedef typename CGAL::Constrained_Delaunay_triangulation_2 CDT; + +typedef CDT::Vertex_handle Vertex_handle; +typedef CDT::Point CDTPoint; + +typedef CGAL::Ray_2 CGAL_Ray_2; +typedef CGAL::Line_3 CGAL_Line_3; +typedef CGAL::Point_2 CGAL_Point_2; +typedef CGAL::Vector_2 CGAL_Vector_2; +typedef CGAL::Segment_2 CGAL_Segment_2; +typedef std::vector CGAL_Polygon_3; +typedef CGAL::Direction_2 CGAL_Direction_2; +typedef CGAL::Direction_3 CGAL_Direction_3; + +/* The idea of 'projection' is how we make 3d points appear as though +they were 2d points to the tessellation algorithm. We take the 3-d plane +on which the polygon lies, and then 'project' or 'cast its shadow' onto +one of three standard planes, the xyplane, the yzplane, or the xzplane, +depending on which projection will prevent the polygon looking like a +flat line. (imagine, the triangle 0,0,1 0,1,1 0,1,0 ... if viewed from +the 'top' it looks line a flat line. so we want to view it from the +side). Thus we create a sequence of x,y points to feed to the algorithm, +but those points might actually be x,z pairs or y,z pairs... it is an +illusion we present to the triangulation algorithm by way of 'projection'. +We get a resulting sequence of triangles with x,y coordinates, which we +then 'deproject' back to x,z or y,z, in 3d space as needed. For example +the square 0,0,0 0,0,1 0,1,1 0,1,0 becomes '0,0 0,1 1,1 1,0', is then +split into two triangles, 0,0 1,0 1,1 and 0,0 1,1 0,1. those two triangles +then are projected back to 3d as 0,0,0 0,1,0 0,1,1 and 0,0 0,1,1 0,0,1. + +There is an additional trick we do with projection related to Polygon +orientation and the orientation of our output triangles, and thus, which +way they are facing in space (aka their 'normals' or 'oriented side'). + +The basic issues is this: every 3d flat polygon can be thought of as +having two sides. In Computer Graphics the convention is that the +'outside' or 'oriented side' or 'normal' is determined by looking at the +triangle in terms of the 'ordering' or 'winding' of the points. If the +points come in a 'clockwise' order, you must be looking at the triangle +from 'inside'. If the points come in a 'counterclockwise' order, you +must be looking at the triangle from the outside. For example, the +triangle 0,0,0 1,0,0 0,1,0, when viewed from the 'top', has points in a +counterclockwise order, so the 'up' side is the 'normal' or 'outside'. +if you look at that same triangle from the 'bottom' side, the points +will appear to be 'clockwise', so the 'down' side is the 'inside', and is the +opposite of the 'normal' side. + +How do we keep track of all that when doing a triangulation? We could +check each triangle as it was generated, and fix it's orientation before +we feed it back to our output list. That is done by, for example, checking +the orientation of the input polygon and then forcing the triangle to +match that orientation during output. This is what CGAL's Nef Polyhedron +does, you can read it inside /usr/include/CGAL/Nef_polyhedron_3.h. + +Or.... we could actually add an additional 'projection' to the incoming +polygon points so that our triangulation algorithm is guaranteed to +create triangles with the proper orientation in the first place. How? +First, we assume that the triangulation algorithm will always produce +'counterclockwise' triangles in our plain old x-y plane. + +The method is based on the following curious fact: That is, if you take +the points of a polygon, and flip the x,y coordinate of each point, +making y:=x and x:=y, then you essentially get a 'mirror image' of the +original polygon... but the orientation will be flipped. Given a +clockwise polygon, the 'flip' will result in a 'counterclockwise' +polygon mirror-image and vice versa. + +Now, there is a second curious fact that helps us here. In 3d, we are +using the plane equation of ax+by+cz+d=0, where a,b,c determine its +direction. If you notice, there are actually mutiple sets of numbers +a:b:c that will describe the exact same plane. For example the 'ground' +plane, called the XYplane, where z is everywhere 0, has the equation +0x+0y+1z+0=0, simplifying to a solution for x,y,z of z=0 and x,y = any +numbers in your number system. However you can also express this as +0x+0y+-1z=0. The x,y,z solution is the same: z is everywhere 0, x and y +are any number, even though a,b,c are different. We can say that the +plane is 'oriented' differently, if we wish. + +But how can we link that concept to the points on the polygon? Well, if +you generate a plane using the standard plane-equation generation +formula, given three points M,N,P, then you will get a plane equation +with . However if you feed the points in the reverse order, +P,N,M, so that they are now oriented in the opposite order, you will get +a plane equation with the signs flipped. <-a:-b:-c:-d> This means you +can essentially consider that a plane has an 'orientation' based on it's +equation, by looking at the signs of a,b,c relative to some other +quantity. + +This means that you can 'flip' the projection of the input polygon +points so that the projection will match the orientation of the input +plane, thus guaranteeing that the output triangles will be oriented in +the same direction as the input polygon was. In other words, even though +we technically 'lose information' when we project from 3d->2d, we can +actually keep the concept of 'orientation' through the whole +triangulation process, and not have to recalculate the proper +orientation during output. + +For example take two side-squares of a cube and the plane equations +formed by feeding the points in counterclockwise, as if looking in from +outside the cube: + + 0,0,0 0,1,0 0,1,1 0,0,1 <-1:0:0:0> + 1,0,0 1,1,0 1,1,1 1,0,1 <1:0:0:1> + +They are both projected onto the YZ plane. They look the same: + 0,0 1,0 1,1 0,1 + 0,0 1,0 1,1 0,1 + +But the second square plane has opposite orientation, so we flip the x +and y for each point: + 0,0 1,0 1,1 0,1 + 0,0 0,1 1,1 1,0 + +Only now do we feed these two 2-d squares to the tessellation algorithm. +The result is 4 triangles. When de-projected back to 3d, they will have +the appropriate winding that will match that of the original 3d faces. +And the first two triangles will have opposite orientation from the last two. +*/ + +typedef enum { XYPLANE, YZPLANE, XZPLANE, NONE } plane_t; +struct projection_t { + plane_t plane; + bool flip; +}; + +CGAL_Point_2 get_projected_point( CGAL_Point_3 &p3, projection_t projection ) { + NT3 x,y; + if (projection.plane == XYPLANE) { x = p3.x(); y = p3.y(); } + else if (projection.plane == XZPLANE) { x = p3.x(); y = p3.z(); } + else if (projection.plane == YZPLANE) { x = p3.y(); y = p3.z(); } + else if (projection.plane == NONE) { x = 0; y = 0; } + if (projection.flip) return CGAL_Point_2( y,x ); + return CGAL_Point_2( x,y ); +} + +/* given 2d point, 3d plane, and 3d->2d projection, 'deproject' from + 2d back onto a point on the 3d plane. true on failure, false on success */ +bool deproject( CGAL_Point_2 &p2, projection_t &projection, CGAL_Plane_3 &plane, CGAL_Point_3 &p3 ) +{ + NT3 x,y; + CGAL_Line_3 l; + CGAL_Point_3 p; + CGAL_Point_2 pf( p2.x(), p2.y() ); + if (projection.flip) pf = CGAL_Point_2( p2.y(), p2.x() ); + if (projection.plane == XYPLANE) { + p = CGAL_Point_3( pf.x(), pf.y(), 0 ); + l = CGAL_Line_3( p, CGAL_Direction_3(0,0,1) ); + } else if (projection.plane == XZPLANE) { + p = CGAL_Point_3( pf.x(), 0, pf.y() ); + l = CGAL_Line_3( p, CGAL_Direction_3(0,1,0) ); + } else if (projection.plane == YZPLANE) { + p = CGAL_Point_3( 0, pf.x(), pf.y() ); + l = CGAL_Line_3( p, CGAL_Direction_3(1,0,0) ); + } + CGAL::Object obj = CGAL::intersection( l, plane ); + const CGAL_Point_3 *point_test = CGAL::object_cast(&obj); + if (point_test) { + p3 = *point_test; + return false; + } + PRINT("ERROR: deproject failure"); + return true; +} + +/* this simple criteria guarantees CGALs triangulation algorithm will +terminate (i.e. not lock up and freeze the program) */ +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(); } +}; + +NT3 sign( const NT3 &n ) +{ + if (n>0) return NT3(1); + if (n<0) return NT3(-1); + return NT3(0); +} + +/* wedge, also related to 'determinant', 'signed parallelogram area', +'side', 'turn', 'winding', '2d portion of cross-product', etc etc. this +function can tell you whether v1 is 'counterclockwise' or 'clockwise' +from v2, based on the sign of the result. when the input Vectors are +formed from three points, A-B and B-C, it can tell you if the path along +the points A->B->C is turning left or right.*/ +NT3 wedge( CGAL_Vector_2 &v1, CGAL_Vector_2 &v2 ) { + return v1.x()*v2.y()-v2.x()*v1.y(); +} + +/* given a point and a possibly non-simple polygon, determine if the +point is inside the polygon or not, using the given winding rule. note +that even_odd is not implemented. */ +typedef enum { NONZERO_WINDING, EVEN_ODD } winding_rule_t; +bool inside(CGAL_Point_2 &p1,std::vector &pgon, winding_rule_t winding_rule) +{ + NT3 winding_sum = NT3(0); + CGAL_Point_2 p2; + CGAL_Ray_2 eastray(p1,CGAL_Direction_2(1,0)); + for (size_t i=0;i(&obj); + if (point_test) { + p2 = *point_test; + CGAL_Vector_2 v1( p1, p2 ); + CGAL_Vector_2 v2( p2, head ); + NT3 this_winding = wedge( v1, v2 ); + winding_sum += sign(this_winding); + } else { + continue; + } + } + if (winding_sum != NT3(0) && winding_rule == NONZERO_WINDING ) return true; + return false; +} + +projection_t find_good_projection( CGAL_Plane_3 &plane ) +{ + projection_t goodproj; + goodproj.plane = NONE; + goodproj.flip = false; + NT3 qxy = plane.a()*plane.a()+plane.b()*plane.b(); + NT3 qyz = plane.b()*plane.b()+plane.c()*plane.c(); + NT3 qxz = plane.a()*plane.a()+plane.c()*plane.c(); + NT3 min = std::min(qxy,std::min(qyz,qxz)); + if (min==qxy) { + goodproj.plane = XYPLANE; + if (sign(plane.c())>0) goodproj.flip = true; + } else if (min==qyz) { + goodproj.plane = YZPLANE; + if (sign(plane.a())>0) goodproj.flip = true; + } else if (min==qxz) { + goodproj.plane = XZPLANE; + if (sign(plane.b())<0) goodproj.flip = true; + } else PRINT("ERROR: failed to find projection"); + return goodproj; +} + +/* given a single near-planar 3d polygon with holes, tessellate into a +sequence of polygons without holes. as of writing, this means conversion +into a sequence of 3d triangles. the given plane should be the same plane +holding the polygon and it's holes. */ +bool tessellate_3d_face_with_holes( std::vector &polygons, std::vector &triangles, CGAL_Plane_3 &plane ) +{ + if (polygons.size()==1 && polygons[0].size()==3) { + PRINT("input polygon has 3 points. shortcut tessellation."); + CGAL_Polygon_3 t; + t.push_back(polygons[0][2]); + t.push_back(polygons[0][1]); + t.push_back(polygons[0][0]); + triangles.push_back( t ); + return false; + } + bool err = false; + CDT cdt; + std::map vertmap; + + PRINT("finding good projection"); + projection_t goodproj = find_good_projection( plane ); + + PRINTB("plane %s",plane ); + PRINTB("proj: %i %i",goodproj.plane % goodproj.flip); + PRINT("Inserting points and edges into Constrained Delaunay Triangulation"); + std::vector< std::vector > polygons2d; + for (size_t i=0;i vhandles; + std::vector polygon2d; + for (size_t j=0;j list_of_seeds; + for (size_t i=1;i &pgon = polygons2d[i]; + for (size_t j=0;j::iterator li = list_of_seeds.begin(); + for (;li!=list_of_seeds.end();li++) { + //PRINTB("seed %s",*li); + double x = CGAL::to_double( li->x() ); + double y = CGAL::to_double( li->y() ); + PRINTB("seed %f,%f",x%y); + } + PRINT("seeding done"); + + PRINT( "meshing" ); + CGAL::refine_Delaunay_mesh_2_without_edge_refinement( cdt, + list_of_seeds.begin(), list_of_seeds.end(), + DummyCriteria() ); + + PRINT("meshing done"); + // this fails because it calls is_simple and is_simple fails on many + // Nef Polyhedron faces + //CGAL::Orientation original_orientation = + // CGAL::orientation_2( orienpgon.begin(), orienpgon.end() ); + + 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]; + CGAL_Point_3 cp1,cp2,cp3; + CGAL_Polygon_3 pgon; + if (vertmap.count(p1)) cp1 = vertmap[p1]; + else err = deproject( p1, goodproj, plane, cp1 ); + if (vertmap.count(p2)) cp2 = vertmap[p2]; + else err = deproject( p2, goodproj, plane, cp2 ); + if (vertmap.count(p3)) cp3 = vertmap[p3]; + else err = deproject( p3, goodproj, plane, cp3 ); + if (err) PRINT("WARNING: 2d->3d deprojection failure"); + pgon.push_back( cp1 ); + pgon.push_back( cp2 ); + pgon.push_back( cp3 ); + triangles.push_back( pgon ); + } + } + + PRINTB("built %i triangles\n",triangles.size()); + return err; +} +/////// Tessellation end + +/* + Create a PolySet from a Nef Polyhedron 3. return false on success, + true on failure. The trick to this is that Nef Polyhedron3 faces have + 'holes' in them. . . while PolySet (and many other 3d polyhedron + formats) do not allow for holes in their faces. The function documents + the method used to deal with this +*/ +bool createPolySetFromNefPolyhedron3(const CGAL_Nef_polyhedron3 &N, PolySet &ps) +{ + bool err = false; + CGAL_Nef_polyhedron3::Halffacet_const_iterator hfaceti; + CGAL_forall_halffacets( hfaceti, N ) { + CGAL_Plane_3 plane( hfaceti->plane() ); + std::vector polygons; + // the 0-mark-volume is the 'empty' volume of space. skip it. + if (hfaceti->incident_volume()->mark() == 0) continue; + CGAL_Nef_polyhedron3::Halffacet_cycle_const_iterator cyclei; + CGAL_forall_facet_cycles_of( cyclei, hfaceti ) { + CGAL_Nef_polyhedron3::SHalfedge_around_facet_const_circulator c1(cyclei); + CGAL_Nef_polyhedron3::SHalfedge_around_facet_const_circulator c2(c1); + CGAL_Polygon_3 polygon; + CGAL_For_all( c1, c2 ) { + CGAL_Point_3 p = c1->source()->center_vertex()->point(); + polygon.push_back( p ); + } + polygons.push_back( polygon ); + } + + /* at this stage, we have a sequence of polygons. the first + is the "outside edge' or 'body' or 'border', and the rest of the + polygons are 'holes' within the first. there are several + options here to get rid of the holes. we choose to go ahead + and let the tessellater deal with the holes, and then + just output the resulting 3d triangles*/ + std::vector triangles; + bool err = tessellate_3d_face_with_holes( polygons, triangles, plane ); + if (!err) for (size_t i=0;i(root_geom)) { switch (format) { case OPENSCAD_STL: - export_stl(ps, output); + export_stl(*ps, output); break; default: assert(false && "Unsupported file format"); @@ -72,16 +73,17 @@ void exportFile(const class Geometry *root_geom, std::ostream &output, FileForma } } -void export_stl(const PolySet *ps, std::ostream &output) +void export_stl(const PolySet &ps, std::ostream &output) { PolySet triangulated(3); - PolysetUtils::tessellate_faces(*ps, triangulated); + PolysetUtils::tessellate_faces(ps, triangulated); setlocale(LC_NUMERIC, "C"); // Ensure radix is . (not ,) in output output << "solid OpenSCAD_Model\n"; BOOST_FOREACH(const PolySet::Polygon &p, triangulated.polygons) { output << " facet normal 0 0 0\n"; output << " outer loop\n"; + assert(p.size() == 3); // STL only allows triangles BOOST_FOREACH(const Vector3d &v, p) { output << "vertex " << v[0] << " " << v[1] << " " << v[2] << "\n"; } @@ -93,24 +95,11 @@ void export_stl(const PolySet *ps, std::ostream &output) } /*! - Saves the current 3D CGAL Nef polyhedron as STL to the given file. + Saves the given CGAL Polyhedon2 as STL to the given file. The file must be open. */ -void export_stl(const CGAL_Nef_polyhedron *root_N, std::ostream &output) +static void export_stl(const CGAL_Polyhedron &P, std::ostream &output) { - if (!root_N->p3->is_simple()) { - PRINT("Object isn't a valid 2-manifold! Modify your design.\n"); - } - CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); - try { - CGAL_Polyhedron P; - //root_N->p3->convert_to_Polyhedron(P); - bool err = nefworkaround::convert_to_Polyhedron( *(root_N->p3), P ); - if (err) { - PRINT("ERROR: CGAL NefPolyhedron->Polyhedron conversion failed"); - return; - } - typedef CGAL_Polyhedron::Vertex Vertex; typedef CGAL_Polyhedron::Vertex_const_iterator VCI; typedef CGAL_Polyhedron::Facet_const_iterator FCI; @@ -176,15 +165,47 @@ void export_stl(const CGAL_Nef_polyhedron *root_N, std::ostream &output) output << "endsolid OpenSCAD_Model\n"; setlocale(LC_NUMERIC, ""); // Set default locale +} +/*! + Saves the current 3D CGAL Nef polyhedron as STL to the given file. + The file must be open. + */ +void export_stl(const CGAL_Nef_polyhedron *root_N, std::ostream &output) +{ + if (!root_N->p3->is_simple()) { + PRINT("Object isn't a valid 2-manifold! Modify your design.\n"); } - catch (const CGAL::Assertion_exception &e) { - PRINTB("CGAL error in CGAL_Nef_polyhedron3::convert_to_Polyhedron(): %s", e.what()); + + bool usePolySet = false; + if (usePolySet) { + PolySet ps(3); + bool err = createPolySetFromNefPolyhedron3(*(root_N->p3), ps); + if (err) { PRINT("ERROR: Nef->PolySet failed"); } + else { + export_stl(ps, output); + } } - catch (...) { - PRINT("CGAL unknown error in CGAL_Nef_polyhedron3::convert_to_Polyhedron()"); + else { + CGAL::Failure_behaviour old_behaviour = CGAL::set_error_behaviour(CGAL::THROW_EXCEPTION); + try { + CGAL_Polyhedron P; + //root_N->p3->convert_to_Polyhedron(P); + bool err = nefworkaround::convert_to_Polyhedron( *(root_N->p3), P ); + if (err) { + PRINT("ERROR: CGAL NefPolyhedron->Polyhedron conversion failed"); + return; + } + export_stl(P, output); + } + catch (const CGAL::Assertion_exception &e) { + PRINTB("CGAL error in CGAL_Nef_polyhedron3::convert_to_Polyhedron(): %s", e.what()); + } + catch (...) { + PRINT("CGAL unknown error in CGAL_Nef_polyhedron3::convert_to_Polyhedron()"); + } + CGAL::set_error_behaviour(old_behaviour); } - CGAL::set_error_behaviour(old_behaviour); } void export_off(const CGAL_Nef_polyhedron *root_N, std::ostream &output) diff --git a/src/export.h b/src/export.h index 0e75d0a8..28c1c07d 100644 --- a/src/export.h +++ b/src/export.h @@ -17,7 +17,7 @@ void exportFile(const class Geometry *root_geom, std::ostream &output, FileForma void export_png(const class Geometry *root_geom, Camera &c, std::ostream &output); void export_stl(const class CGAL_Nef_polyhedron *root_N, std::ostream &output); -void export_stl(const class PolySet *ps, std::ostream &output); +void export_stl(const class PolySet &ps, std::ostream &output); void export_off(const CGAL_Nef_polyhedron *root_N, std::ostream &output); void export_dxf(const class Polygon2d &poly, std::ostream &output); void export_png(const CGAL_Nef_polyhedron *root_N, Camera &c, std::ostream &output);