diff --git a/openscad.pro b/openscad.pro index f926ac56..6ee6fa34 100644 --- a/openscad.pro +++ b/openscad.pro @@ -440,6 +440,7 @@ HEADERS += src/cgal.h \ SOURCES += src/cgalutils.cc \ src/cgalutils-tess.cc \ + src/cgalutils-tess-old.cc \ src/CGALCache.cc \ src/CGALRenderer.cc \ src/CGAL_Nef_polyhedron.cc \ diff --git a/src/cgalutils-tess-old.cc b/src/cgalutils-tess-old.cc new file mode 100644 index 00000000..d779772f --- /dev/null +++ b/src/cgalutils-tess-old.cc @@ -0,0 +1,553 @@ +/* + +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 "cgalutils.h" +#include +#include + +typedef CGAL_Kernel3 Kernel; +//typedef CGAL::Triangulation_vertex_base_2 Vb; +typedef CGAL::Triangulation_vertex_base_2 Vb; +//typedef CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Delaunay_mesh_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 TDS; +typedef CGAL::Exact_intersections_tag ITAG; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +//typedef 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 CGAL::Direction_2 CGAL_Direction_2; +typedef CGAL::Direction_3 CGAL_Direction_3; +typedef CGAL::Plane_3 CGAL_Plane_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; +} + +namespace CGALUtils { +/* 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 tessellate3DFaceWithHolesNew(std::vector &polygons, + Polygons &triangles, + CGAL_Plane_3 &plane) + { + if (polygons.size()==1 && polygons[0].size()==3) { + PRINTD("input polygon has 3 points. shortcut tessellation."); + Polygon t; + t.push_back(Vector3d(CGAL::to_double(polygons[0][0].x()), CGAL::to_double(polygons[0][0].y()), CGAL::to_double(polygons[0][0].z()))); + t.push_back(Vector3d(CGAL::to_double(polygons[0][1].x()), CGAL::to_double(polygons[0][1].y()), CGAL::to_double(polygons[0][1].z()))); + t.push_back(Vector3d(CGAL::to_double(polygons[0][2].x()), CGAL::to_double(polygons[0][2].y()), CGAL::to_double(polygons[0][2].z()))); + triangles.push_back( t ); + return false; + } + bool err = false; + CDT cdt; + std::map vertmap; + + PRINTD("finding good projection"); + projection_t goodproj = find_good_projection( plane ); + + PRINTDB("plane %s",plane ); + PRINTDB("proj: %i %i",goodproj.plane % goodproj.flip); + PRINTD("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() ); + PRINTDB("seed %f,%f",x%y); + } + PRINTD("seeding done"); + + PRINTD( "meshing" ); + CGAL::refine_Delaunay_mesh_2_without_edge_refinement( cdt, + list_of_seeds.begin(), list_of_seeds.end(), + DummyCriteria() ); + + PRINTD("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; + 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"); + Polygon tri; + tri.push_back(Vector3d(CGAL::to_double(cp1.x()), CGAL::to_double(cp1.y()), CGAL::to_double(cp1.z()))); + tri.push_back(Vector3d(CGAL::to_double(cp2.x()), CGAL::to_double(cp2.y()), CGAL::to_double(cp2.z()))); + tri.push_back(Vector3d(CGAL::to_double(cp3.x()), CGAL::to_double(cp3.y()), CGAL::to_double(cp3.z()))); + triangles.push_back( tri ); + } + } + + PRINTDB("built %i triangles",triangles.size()); + return err; + } + + +/* 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 tessellate3DFaceWithHoles(std::vector &polygons, + std::vector &triangles, + CGAL_Plane_3 &plane) + { + if (polygons.size()==1 && polygons[0].size()==3) { + PRINTD("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; + + PRINTD("finding good projection"); + projection_t goodproj = find_good_projection( plane ); + + PRINTDB("plane %s",plane ); + PRINTDB("proj: %i %i",goodproj.plane % goodproj.flip); + PRINTD("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() ); + PRINTDB("seed %f,%f",x%y); + } + PRINTD("seeding done"); + + PRINTD( "meshing" ); + CGAL::refine_Delaunay_mesh_2_without_edge_refinement( cdt, + list_of_seeds.begin(), list_of_seeds.end(), + DummyCriteria() ); + + PRINTD("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 ); + } + } + + PRINTDB("built %i triangles",triangles.size()); + return err; + } + +}; diff --git a/src/cgalutils-tess.cc b/src/cgalutils-tess.cc index c8188908..e230e9db 100644 --- a/src/cgalutils-tess.cc +++ b/src/cgalutils-tess.cc @@ -1,433 +1,198 @@ -/* - -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 "cgalutils.h" -#include -#include +//#include "cgal.h" +//#include "tess.h" -typedef CGAL_Kernel3 Kernel; -//typedef CGAL::Triangulation_vertex_base_2 Vb; -typedef CGAL::Triangulation_vertex_base_2 Vb; -//typedef CGAL::Constrained_triangulation_face_base_2 Fb; -typedef CGAL::Delaunay_mesh_face_base_2 Fb; -typedef CGAL::Triangulation_data_structure_2 TDS; -typedef CGAL::Exact_intersections_tag ITAG; -typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; -//typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +#ifdef NDEBUG +#define PREV_NDEBUG NDEBUG +#undef NDEBUG +#endif +#include +#include +#include +#ifdef PREV_NDEBUG +#define NDEBUG PREV_NDEBUG +#endif -typedef CDT::Vertex_handle Vertex_handle; -typedef CDT::Point CDTPoint; +#include -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 CGAL::Direction_2 CGAL_Direction_2; -typedef CGAL::Direction_3 CGAL_Direction_3; -typedef CGAL::Plane_3 CGAL_Plane_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; +struct FaceInfo { + int nesting_level; + bool in_domain() { return nesting_level%2 == 1; } }; -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 ); -} +typedef CGAL::Triangulation_2_filtered_projection_traits_3 Projection; +typedef CGAL::Triangulation_face_base_with_info_2 Fbb; +typedef CGAL::Triangulation_data_structure_2< + CGAL::Triangulation_vertex_base_2, + CGAL::Constrained_triangulation_face_base_2 > Tds; +typedef CGAL::Constrained_Delaunay_triangulation_2< + Projection, Tds, CGAL::Exact_predicates_tag> CDT; -/* 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 ) + +static void mark_domains(CDT &ct, + CDT::Face_handle start, + int index, + std::list& border) { - 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; + if (start->info().nesting_level != -1) return; + std::list queue; + queue.push_back(start); + while (!queue.empty()) { + CDT::Face_handle fh = queue.front(); + queue.pop_front(); + if (fh->info().nesting_level == -1) { + fh->info().nesting_level = index; + for (int i = 0; i < 3; i++) { + CDT::Edge e(fh,i); + CDT::Face_handle n = fh->neighbor(i); + if (n->info().nesting_level == -1) { + if (ct.is_constrained(e)) border.push_back(e); + else queue.push_back(n); + } + } + } + } } -/* 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 ) +//explore set of facets connected with non constrained edges, +//and attribute to each such set a nesting level. +//We start from facets incident to the infinite vertex, with a nesting +//level of 0. Then we recursively consider the non-explored facets incident +//to constrained edges bounding the former set and increase the nesting level by 1. +//Facets in the domain are those with an odd nesting level. +static void mark_domains(CDT& cdt) { - 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; + for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it) { + it->info().nesting_level = -1; + } + std::list border; + mark_domains(cdt, cdt.infinite_face(), 0, border); + while (!border.empty()) { + CDT::Edge e = border.front(); + border.pop_front(); + CDT::Face_handle n = e.first->neighbor(e.second); + if (n->info().nesting_level == -1) { + mark_domains(cdt, n, e.first->info().nesting_level+1, border); + } + } } namespace CGALUtils { -/* 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 tessellate3DFaceWithHoles(std::vector &polygons, - std::vector &triangles, - CGAL_Plane_3 &plane) + /*! + polygons define a polygon, potentially with holes. Each contour + should be added as a separate PolygonK instance. + The tessellator will handle almost planar polygons. + + If the normal is given, we will assume this as the normal vector of the polygon. + Otherwise, we will try to calculate it using Newell's method. + + The resulting triangles is added to the given triangles vector. + */ + bool tessellatePolygonWithHoles(const PolyholeK &polygons, + Polygons &triangles, + const K::Vector_3 *normal) { - if (polygons.size()==1 && polygons[0].size()==3) { + // No polygon. FIXME: Will this ever happen or can we assert here? + if (polygons.empty()) return false; + + // No hole + if (polygons.size() == 1) return tessellatePolygon(polygons.front(), triangles, normal); + + K::Vector_3 normalvec; + if (normal) { + normalvec = *normal; + } + else { + // Calculate best guess at face normal using Newell's method + CGAL::normal_vector_newell_3(polygons.front().begin(), polygons.front().end(), normalvec); + } + + // Pass the normal vector to the (undocumented) + // CGAL::Triangulation_2_filtered_projection_traits_3. This + // trait deals with projection from 3D to 2D using the normal + // vector as a hint, and allows for near-planar polygons to be passed to + // the Constrained Delaunay Triangulator. + Projection actualProjection(normalvec); + CDT cdt(actualProjection); + BOOST_FOREACH(const PolygonK &poly, polygons) { + for (size_t i=0;iinfo().in_domain()) { + Polygon tri; + for (int i=0;i<3;i++) { + Vertex3K v = cdt.triangle(fit)[i]; + tri.push_back(Vector3d(v.x(), v.y(), v.z())); + } + triangles.push_back(tri); + } + } + + return false; + } + + bool tessellatePolygon(const PolygonK &polygon, + Polygons &triangles, + const K::Vector_3 *normal) + { + if (polygon.size() == 3) { PRINTD("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 ); + Polygon t; + t.push_back(Vector3d(polygon[0].x(), polygon[0].y(), polygon[0].z())); + t.push_back(Vector3d(polygon[1].x(), polygon[1].y(), polygon[1].z())); + t.push_back(Vector3d(polygon[2].x(), polygon[2].y(), polygon[2].z())); + triangles.push_back(t); return false; } - bool err = false; - CDT cdt; - std::map vertmap; - PRINTD("finding good projection"); - projection_t goodproj = find_good_projection( plane ); - - PRINTDB("plane %s",plane ); - PRINTDB("proj: %i %i",goodproj.plane % goodproj.flip); - PRINTD("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() ); - PRINTDB("seed %f,%f",x%y); + + // Pass the normal vector to the (undocumented) + // CGAL::Triangulation_2_filtered_projection_traits_3. This + // trait deals with projection from 3D to 2D using the normal + // vector as a hint, and allows for near-planar polygons to be passed to + // the Constrained Delaunay Triangulator. + Projection actualProjection(normalvec); + CDT cdt(actualProjection); + for (size_t i=0;i() ); - - PRINTD("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() ); + + //Mark facets that are inside the domain bounded by the polygon + mark_domains(cdt); + // Iterate over the resulting faces 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 ); + for (fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); fit++) { + if (fit->info().in_domain()) { + Polygon tri; + for (int i=0;i<3;i++) { + K::Point_3 v = cdt.triangle(fit)[i]; + tri.push_back(Vector3d(v.x(), v.y(), v.z())); + } + triangles.push_back(tri); } } - PRINTDB("built %i triangles",triangles.size()); - return err; + return false; } -}; + +}; // namespace CGALUtils + diff --git a/src/cgalutils.cc b/src/cgalutils.cc index a609bd63..5338a14e 100644 --- a/src/cgalutils.cc +++ b/src/cgalutils.cc @@ -21,21 +21,21 @@ #include #include +namespace Eigen { + size_t hash_value(Vector3d const &v) { + size_t seed = 0; + for (int i=0;i<3;i++) boost::hash_combine(seed, v[i]); + return seed; + } +} + 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])); } -#undef GEN_SURFACE_DEBUG - - namespace Eigen { - size_t hash_value(Vector3d const &v) { - size_t seed = 0; - for (int i=0;i<3;i++) boost::hash_combine(seed, v[i]); - return seed; - } - } +#define GEN_SURFACE_DEBUG class CGAL_Build_PolySet : public CGAL::Modifier_base { @@ -138,7 +138,6 @@ namespace /* anonymous */ { void operator()(CGAL_HDS& hds) { CGAL_Polybuilder B(hds, true); - typedef boost::tuple BuilderVertex; Reindexer vertices; std::vector indices(3); @@ -230,14 +229,14 @@ namespace /* anonymous */ { 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())); + 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()); + Index index(in_poly.vertices_begin(), in_poly.vertices_end()); for(Facet_const_iterator fi = in_poly.facets_begin(), end = in_poly.facets_end(); @@ -245,13 +244,13 @@ namespace /* anonymous */ { { HFCC hc = fi->facet_begin(); HFCC hc_end = hc; - // std::size_t n = circulator_size( hc); - // CGAL_assertion( n >= 3); + // 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); + } while(hc != hc_end); builder.end_facet(); } builder.end_surface(); @@ -300,19 +299,18 @@ static CGAL_Nef_polyhedron *createNefPolyhedronFromPolySet(const PolySet &ps) try { CGAL_Polyhedron P; bool err = CGALUtils::createPolyhedronFromPolySet(ps, P); - // if (!err) { - // PRINTB("Polyhedron is closed: %d", P.is_closed()); - // PRINTB("Polyhedron is valid: %d", P.is_valid(false, 0)); - // } + if (!err) { + PRINTDB("Polyhedron is closed: %d", P.is_closed()); + PRINTDB("Polyhedron is valid: %d", P.is_valid(false, 0)); + } 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) { + if (std::string(e.what()).find("Plane_constructor")!=std::string::npos && + 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()); } @@ -461,7 +459,7 @@ namespace CGALUtils { if ((ps && ps->is_convex()) || (!ps && is_weakly_convex(poly))) { - PRINTDB("Minkowski: child %d is convex and %s",i % (ps?"PolySet":"Nef") ); + PRINTDB("Minkowski: child %d is convex and %s",i % (ps?"PolySet":"Nef")); P[i].push_back(poly); } else { CGAL_Nef_polyhedron3 decomposed_nef; @@ -480,7 +478,7 @@ namespace CGALUtils { // the first volume is the outer volume, which ignored in the decomposition CGAL_Nef_polyhedron3::Volume_const_iterator ci = ++decomposed_nef.volumes_begin(); - for( ; ci != decomposed_nef.volumes_end(); ++ci) { + for(; ci != decomposed_nef.volumes_end(); ++ci) { if(ci->mark()) { CGAL_Polyhedron poly; decomposed_nef.convert_inner_shell_to_polyhedron(ci->shells_begin(), poly); @@ -844,13 +842,13 @@ namespace CGALUtils { // dont use z of 0. there are bugs in CGAL. double inf = 1e8; double eps = 0.001; - CGAL_Point_3 minpt( -inf, -inf, -eps ); - CGAL_Point_3 maxpt( inf, inf, eps ); - CGAL_Iso_cuboid_3 bigcuboid( minpt, maxpt ); - for ( int i=0;i<8;i++ ) pts.push_back( bigcuboid.vertex(i) ); + CGAL_Point_3 minpt(-inf, -inf, -eps); + CGAL_Point_3 maxpt( inf, inf, eps); + CGAL_Iso_cuboid_3 bigcuboid(minpt, maxpt); + for (int i=0;i<8;i++) pts.push_back(bigcuboid.vertex(i)); CGAL_Polyhedron bigbox; CGAL::convex_hull_3(pts.begin(), pts.end(), bigbox); - CGAL_Nef_polyhedron3 nef_bigbox( bigbox ); + CGAL_Nef_polyhedron3 nef_bigbox(bigbox); newN.p3.reset(new CGAL_Nef_polyhedron3(nef_bigbox.intersection(*N.p3))); } catch (const CGAL::Failure_exception &e) { @@ -864,18 +862,18 @@ namespace CGALUtils { return poly; } - PRINTDB("%s",OpenSCAD::svg_header( 480, 100000 )); + PRINTDB("%s",OpenSCAD::svg_header(480, 100000)); try { ZRemover zremover; CGAL_Nef_polyhedron3::Volume_const_iterator i; CGAL_Nef_polyhedron3::Shell_entry_const_iterator j; CGAL_Nef_polyhedron3::SFace_const_handle sface_handle; - for ( i = newN.p3->volumes_begin(); i != newN.p3->volumes_end(); ++i ) { + for (i = newN.p3->volumes_begin(); i != newN.p3->volumes_end(); ++i) { PRINTDB("",i->mark()); - for ( j = i->shells_begin(); j != i->shells_end(); ++j ) { + for (j = i->shells_begin(); j != i->shells_end(); ++j) { PRINTDB(""); } PRINTD(""); @@ -909,7 +907,7 @@ namespace CGALUtils { // can be optimized by rewriting bounding_box to accept vertices CGAL_forall_vertices(vi, N) points.push_back(vi->point()); - if (points.size()) result = CGAL::bounding_box( points.begin(), points.end() ); + if (points.size()) result = CGAL::bounding_box(points.begin(), points.end()); return result; } @@ -1052,25 +1050,26 @@ namespace CGALUtils { formats) do not allow for holes in their faces. The function documents the method used to deal with this */ +#if 0 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; + CGAL_forall_halffacets(hfaceti, N) { + CGAL::Plane_3 plane(hfaceti->plane()); + std::vector polyholes; // the 0-mark-volume is the 'empty' volume of space. skip it. if (hfaceti->incident_volume()->mark()) continue; CGAL_Nef_polyhedron3::Halffacet_cycle_const_iterator cyclei; - CGAL_forall_facet_cycles_of( cyclei, hfaceti ) { + 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_For_all(c1, c2) { CGAL_Point_3 p = c1->source()->center_vertex()->point(); - polygon.push_back( p ); + polygon.push_back(p); } - polygons.push_back( polygon ); + polyholes.push_back(polygon); } /* at this stage, we have a sequence of polygons. the first @@ -1079,28 +1078,125 @@ namespace CGALUtils { 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*/ + + CGAL::Vector_3 nvec = plane.orthogonal_vector(); + K::Vector_3 normal(CGAL::to_double(nvec.x()), CGAL::to_double(nvec.y()), CGAL::to_double(nvec.z())); std::vector triangles; - bool err = CGALUtils::tessellate3DFaceWithHoles(polygons, triangles, plane); + bool err = CGALUtils::tessellate3DFaceWithHoles(polyholes, triangles, plane); if (!err) { - for (size_t i=0;i=0;j--) { - double x1,y1,z1; - x1 = CGAL::to_double(triangles[i][j].x()); - y1 = CGAL::to_double(triangles[i][j].y()); - z1 = CGAL::to_double(triangles[i][j].z()); - ps.append_vertex(x1,y1,z1); - } + ps.append_vertex(CGAL::to_double(p[2].x()), CGAL::to_double(p[2].y()), CGAL::to_double(p[2].z())); + ps.append_vertex(CGAL::to_double(p[1].x()), CGAL::to_double(p[1].y()), CGAL::to_double(p[1].z())); + ps.append_vertex(CGAL::to_double(p[0].x()), CGAL::to_double(p[0].y()), CGAL::to_double(p[0].z())); } } } return err; } +#endif +#if 0 + 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 polyholes; + // the 0-mark-volume is the 'empty' volume of space. skip it. + if (hfaceti->incident_volume()->mark()) 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); + } + polyholes.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*/ + + CGAL::Vector_3 nvec = plane.orthogonal_vector(); + K::Vector_3 normal(CGAL::to_double(nvec.x()), CGAL::to_double(nvec.y()), CGAL::to_double(nvec.z())); + std::vector triangles; + bool err = CGALUtils::tessellate3DFaceWithHolesNew(polyholes, triangles, plane); + if (!err) { + BOOST_FOREACH(const Polygon &p, triangles) { + if (p.size() != 3) { + PRINT("WARNING: triangle doesn't have 3 points. skipping"); + continue; + } + ps.append_poly(); + ps.append_vertex(p[0].x(), p[0].y(), p[0].z()); + ps.append_vertex(p[1].x(), p[1].y(), p[1].z()); + ps.append_vertex(p[2].x(), p[2].y(), p[2].z()); + } + } + } + return err; + } +#endif +#if 1 + 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()); + PolyholeK polyholes; + // the 0-mark-volume is the 'empty' volume of space. skip it. + if (hfaceti->incident_volume()->mark()) 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); + PolygonK polygon; + CGAL_For_all(c1, c2) { + CGAL_Point_3 p = c1->source()->center_vertex()->point(); + polygon.push_back(Vertex3K(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z()))); + } + polyholes.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*/ + + CGAL::Vector_3 nvec = plane.orthogonal_vector(); + K::Vector_3 normal(CGAL::to_double(nvec.x()), CGAL::to_double(nvec.y()), CGAL::to_double(nvec.z())); + std::vector triangles; + bool err = CGALUtils::tessellatePolygonWithHoles(polyholes, triangles, &normal); + if (!err) { + BOOST_FOREACH(const Polygon &p, triangles) { + if (p.size() != 3) { + PRINT("WARNING: triangle doesn't have 3 points. skipping"); + continue; + } + ps.append_poly(); + ps.append_vertex(p[0].x(), p[0].y(), p[0].z()); + ps.append_vertex(p[1].x(), p[1].y(), p[1].z()); + ps.append_vertex(p[2].x(), p[2].y(), p[2].z()); + } + } + } + return err; + } +#endif CGAL_Nef_polyhedron *createNefPolyhedronFromGeometry(const Geometry &geom) { const PolySet *ps = dynamic_cast(&geom); @@ -1117,10 +1213,10 @@ namespace CGALUtils { }; // namespace CGALUtils -void ZRemover::visit( CGAL_Nef_polyhedron3::Halffacet_const_handle hfacet ) +void ZRemover::visit(CGAL_Nef_polyhedron3::Halffacet_const_handle hfacet) { PRINTDB(" ",hfacet->mark()); - if ( hfacet->plane().orthogonal_direction() != this->up ) { + if (hfacet->plane().orthogonal_direction() != this->up) { PRINTD(" "); PRINTD(" "); return; @@ -1130,36 +1226,36 @@ void ZRemover::visit( CGAL_Nef_polyhedron3::Halffacet_const_handle hfacet ) CGAL_Nef_polyhedron3::Halffacet_cycle_const_iterator fci; int contour_counter = 0; - CGAL_forall_facet_cycles_of( fci, hfacet ) { - if ( fci.is_shalfedge() ) { + CGAL_forall_facet_cycles_of(fci, hfacet) { + if (fci.is_shalfedge()) { PRINTD(" "); CGAL_Nef_polyhedron3::SHalfedge_around_facet_const_circulator c1(fci), cend(c1); std::vector contour; - CGAL_For_all( c1, cend ) { + CGAL_For_all(c1, cend) { CGAL_Nef_polyhedron3::Point_3 point3d = c1->source()->target()->point(); CGAL_Nef_polyhedron2::Explorer::Point point2d(CGAL::to_double(point3d.x()), CGAL::to_double(point3d.y())); - contour.push_back( point2d ); + contour.push_back(point2d); } if (contour.size()==0) continue; if (OpenSCAD::debug!="") - PRINTDB(" ", CGAL::is_simple_2( contour.begin(), contour.end() ) ); + PRINTDB(" ", CGAL::is_simple_2(contour.begin(), contour.end())); - tmpnef2d.reset( new CGAL_Nef_polyhedron2( contour.begin(), contour.end(), boundary ) ); + tmpnef2d.reset(new CGAL_Nef_polyhedron2(contour.begin(), contour.end(), boundary)); - if ( contour_counter == 0 ) { - PRINTDB(" ", contour.size() ); + if (contour_counter == 0) { + PRINTDB(" ", contour.size()); *(output_nefpoly2d) += *(tmpnef2d); } else { - PRINTDB(" ", contour.size() ); + PRINTDB(" ", contour.size()); *(output_nefpoly2d) *= *(tmpnef2d); } /*log << "\n\n" - << OpenSCAD::dump_svg( *tmpnef2d ) << "\n" + << OpenSCAD::dump_svg(*tmpnef2d) << "\n" << "\n\n" - << OpenSCAD::dump_svg( *output_nefpoly2d ) << "\n";*/ + << OpenSCAD::dump_svg(*output_nefpoly2d) << "\n";*/ contour_counter++; } else { diff --git a/src/cgalutils.h b/src/cgalutils.h index f660f5aa..225e93f1 100644 --- a/src/cgalutils.h +++ b/src/cgalutils.h @@ -5,6 +5,12 @@ #include "CGAL_Nef_polyhedron.h" #include "enums.h" +#include +typedef CGAL::Epick K; +typedef CGAL::Point_3 Vertex3K; +typedef std::vector PolygonK; +typedef std::vector PolyholeK; + namespace CGALUtils { bool applyHull(const Geometry::ChildList &children, PolySet &P); void applyOperator(const Geometry::ChildList &children, CGAL_Nef_polyhedron &dest, OpenSCADOperator op); @@ -19,8 +25,16 @@ namespace CGALUtils { bool createPolySetFromNefPolyhedron3(const CGAL_Nef_polyhedron3 &N, PolySet &ps); bool createPolyhedronFromPolySet(const PolySet &ps, CGAL_Polyhedron &p); - typedef std::vector CGAL_Polygon_3; - bool tessellate3DFaceWithHoles(std::vector &polygons, + bool tessellatePolygon(const PolygonK &polygon, + Polygons &triangles, + const K::Vector_3 *normal = NULL); + bool tessellatePolygonWithHoles(const PolyholeK &polygons, + Polygons &triangles, + const K::Vector_3 *normal = NULL); + bool tessellate3DFaceWithHolesNew(std::vector &polygons, + Polygons &triangles, + CGAL::Plane_3 &plane); + bool tessellate3DFaceWithHoles(std::vector &polygons, std::vector &triangles, CGAL::Plane_3 &plane); }; diff --git a/src/polyset-utils.cc b/src/polyset-utils.cc index b55dfc4d..5587be6e 100644 --- a/src/polyset-utils.cc +++ b/src/polyset-utils.cc @@ -2,25 +2,7 @@ #include "polyset.h" #include "Polygon2d.h" #include "printutils.h" -#include "cgal.h" - -#ifdef NDEBUG -#define PREV_NDEBUG NDEBUG -#undef NDEBUG -#endif -#include -#include -#include -#ifdef PREV_NDEBUG -#define NDEBUG PREV_NDEBUG -#endif - -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Triangulation_2_filtered_projection_traits_3 Projection; -typedef CGAL::Triangulation_data_structure_2 < - CGAL::Triangulation_vertex_base_2, - CGAL::Constrained_triangulation_face_base_2 > Tds; -typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +#include "cgalutils.h" #include @@ -60,7 +42,7 @@ namespace PolysetUtils { 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. - The tessellation will be robust wrt. degenerate and self-intersecting + The tessellation will be robust wrt. degenerate and self-intersecting */ void tessellate_faces(const PolySet &inps, PolySet &outps) { int degeneratePolygons = 0; @@ -76,35 +58,12 @@ namespace PolysetUtils { } else { // Build a data structure that CGAL accepts - std::vector cgalpoints; + PolygonK cgalpoints; BOOST_FOREACH(const Vector3d &v, pgon) { - cgalpoints.push_back(K::Point_3(v[0], v[1], v[2])); + cgalpoints.push_back(Vertex3K(v[0], v[1], v[2])); } - // Calculate best guess at face normal using Newell's method - K::Vector_3 normal; - CGAL::normal_vector_newell_3(cgalpoints.begin(), cgalpoints.end(), normal); - // Pass the normal vector to the (undocumented) - // CGAL::Triangulation_2_filtered_projection_traits_3. This - // trait deals with projection from 3D to 2D using the normal - // vector as a hint, and allows for near-planar polygons to be passed to - // the Constrained Delaunay Triangulator. - Projection actualProjection(normal); - CDT cdt(actualProjection); - for (size_t i=0;i