diff --git a/Packages/Triangulation_3/include/CGAL/Delaunay_remove_tds_3.h b/Packages/Triangulation_3/include/CGAL/Delaunay_remove_tds_3.h index 4c3e33f4f19..0bb336fc23a 100644 --- a/Packages/Triangulation_3/include/CGAL/Delaunay_remove_tds_3.h +++ b/Packages/Triangulation_3/include/CGAL/Delaunay_remove_tds_3.h @@ -31,26 +31,32 @@ public : : _f(f) {} - void* face() const { return _f; + void* + face() const { return _f; } - void set_face(void* f) { + void + set_face(void* f) { _f = f ; } - bool is_valid(bool /* verbose */ = false, int /* level */ = 0) const { + bool + is_valid(bool /* verbose */ = false, int /* level */ = 0) const { return true; } - void set_info(I i) { + void + set_info(I i) { _info = i; } - I info(){ + I + info(){ return _info; } - Point point() { + Point + point() { return _info->point(); } @@ -60,7 +66,15 @@ private: }; -// We derive the face class. +/* We derive the face class, + - because we want an additional pointer, + - because we want to manage the order of the list of the faces + in the triangulation datastructure, + - because we want to look at each 'edge' only once. If two + triangles are neighbors, then the triangle with the smaller + address has the edge. +*/ + template < class Gt,class I > class Delaunay_remove_tds_face_3_2 :public Triangulation_face_base_2 { @@ -75,7 +89,8 @@ public: void* n0, void* n1, void* n2) : Triangulation_face_base_2(v0, v1, v2, n0, n1, n2) {} - void set_info(I i) { + void + set_info(I i) { inf = i; } @@ -84,13 +99,21 @@ public: } //Handling the doubly connected list of faces - void* p() const {return _p;} - void* n() const {return _n;} - void set_p(void* f) {_p = f;}; - void set_n(void* f) {_n = f;}; + void* + p() const {return _p;} + + void* + n() const {return _n;} + + void + set_p(void* f) {_p = f;}; + + void + set_n(void* f) {_n = f;}; // Remove this face from the list - void unlink() { + void + unlink() { if(_p) { ((Delaunay_remove_tds_face_3_2*)_p)->set_n(_n); } @@ -100,7 +123,8 @@ public: } // Remove neighbor i from the list - void unlink(int i) { + void + unlink(int i) { Delaunay_remove_tds_face_3_2 * n = (Delaunay_remove_tds_face_3_2*)neighbor(cw(i)); n->unlink(); @@ -109,7 +133,9 @@ public: } // Move a given face after this - void move_after_this(Delaunay_remove_tds_face_3_2* f) { +private: + void + move_after_this(Delaunay_remove_tds_face_3_2* f) { if(_n == f) { return; } @@ -125,7 +151,11 @@ public: f->set_p(this); } - void set_edge(int i, Delaunay_remove_tds_face_3_2* f) { +public: + // Note that when we set an edge, the face gets moved + // so that it gets considered later. + void + set_edge(int i, Delaunay_remove_tds_face_3_2* f) { Delaunay_remove_tds_face_3_2 *n = (Delaunay_remove_tds_face_3_2*)neighbor(i); if(n < this) { @@ -138,22 +168,26 @@ public: } } - void set_edge(int i) { + void + set_edge(int i) { set_edge(i, this); } - void set_edge() { + void + set_edge() { for(int i = 0; i < 3; i++) { set_edge(i, this); } } - void set_edge(int i, bool b) { + void + set_edge(int i, bool b) { _edge[i] = b; } - bool edge(int i) const { + bool + edge(int i) const { return _edge[i]; } @@ -181,6 +215,7 @@ public: // This class is used to represent the boundary of a hole in a polyhedron. +// It only implements a constructor, the rest is inherited template class Delaunay_remove_tds_3_2 @@ -246,40 +281,29 @@ public: Vertex_3_2 *v0, *v1, *v2; Vertex *w0 = handle2pointer(h->vertex(i0)); - map_it = vertex_map.lower_bound(w0); - if((map_it == vertex_map.end()) || (w0 != map_it->first)) { + + if((map_it = vertex_map.find(w0)) == vertex_map.end()) { v0 = create_vertex(); v0->set_info(w0); - if((map_it != vertex_map.begin()) && (map_it != vertex_map.end())) { - map_it--; - } - vertex_map.insert(map_it, std::make_pair(w0, v0)); + vertex_map.insert(std::make_pair(w0, v0)); } else { v0 = map_it->second; } - + Vertex *w1 = handle2pointer(h->vertex(i1)); - map_it = vertex_map.lower_bound(w1); - if((map_it == vertex_map.end()) || (w1 != map_it->first)) { + if((map_it = vertex_map.find(w1)) == vertex_map.end()) { v1 = create_vertex(); v1->set_info(w1); - if((map_it != vertex_map.begin()) && (map_it != vertex_map.end())) { - map_it--; - } - vertex_map.insert(map_it, std::make_pair(w1, v1)); + vertex_map.insert(std::make_pair(w1, v1)); } else { v1 = map_it->second; } Vertex *w2 = handle2pointer(h->vertex(i2)); - map_it = vertex_map.lower_bound(w2); - if((map_it == vertex_map.end()) || (w2 != map_it->first)) { + if((map_it = vertex_map.find(w2)) == vertex_map.end()) { v2 = create_vertex(); v2->set_info(w2); - if((map_it != vertex_map.begin()) && (map_it != vertex_map.end())) { - map_it--; - } - vertex_map.insert(map_it, std::make_pair(w2, v2)); + vertex_map.insert(std::make_pair(w2, v2)); } else { v2 = map_it->second; } diff --git a/Packages/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h b/Packages/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h index b4b2b35cacb..b699dc6b8fb 100644 --- a/Packages/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h +++ b/Packages/Triangulation_3/include/CGAL/Delaunay_triangulation_3.h @@ -1637,7 +1637,6 @@ remove_3D_ear(Vertex_handle v) CGAL_triangulation_precondition( dimension() == 3 ); std::list boundary; // facets on the boundary of the hole - // std::set bdvert; // vertices on the boundary of the hole if ( test_dim_down(v) ) { @@ -1663,7 +1662,6 @@ remove_3D_ear(Vertex_handle v) std::list hole; make_hole_3D_ear(v, boundary, - //bdvert, hole); @@ -1687,7 +1685,6 @@ void Delaunay_triangulation_3:: make_hole_3D_ear( Vertex_handle v, std::list & boundhole, - //std::set & boundvert, std::list & hole) { CGAL_triangulation_precondition( ! test_dim_down(v) ); @@ -1712,11 +1709,7 @@ make_hole_3D_ear( Vertex_handle v, for ( i=0; i<4; i++) { if ( i != indv ) { vi = (*cit)->vertex(i); - //if ( boundvert.find( vi ) == boundvert.end() ) - // { - // boundvert.insert( vi ); vi->set_cell( opp_cit ); - // } } } @@ -1731,7 +1724,7 @@ template < class Gt, class Tds > void Delaunay_triangulation_3:: undo_make_hole_3D_ear(std::list & boundhole, - std::list & hole){ + std::list & hole){ std::list::iterator cit = hole.begin(); for(std::list::iterator fit = boundhole.begin(); fit != boundhole.end(); @@ -1767,7 +1760,7 @@ fill_hole_3D_ear( std::list & boundhole) int opcount = 0; Face_3_2 *f = &(* surface.faces_begin()); - Face_3_2 *last_op; + Face_3_2 *last_op = NULL; // This is where the last ear was inserted int k = -1; // This is a loop over the halfedges of the surface of the hole @@ -1783,8 +1776,8 @@ fill_hole_3D_ear( std::list & boundhole) } k = 0; } + if(f->edge(k)) { - bool violation = false; Vertex_3_2 *w0, *w1, *w2, *w3; Vertex *v0, *v1, *v2, *v3; int i = ccw(k); @@ -1807,7 +1800,6 @@ fill_hole_3D_ear( std::list & boundhole) const Point & p2 = v2->point(); const Point & p3 = v3->point(); - bool inf_0 = is_infinite(Vertex_handle(v0)); bool inf_1 = is_infinite(Vertex_handle(v1)); bool inf_2 = is_infinite(Vertex_handle(v2)); @@ -1816,8 +1808,7 @@ fill_hole_3D_ear( std::list & boundhole) if( inf_1 || inf_2 ){ // there will be another ear, so let's ignore this one, // because it is complicated to treat - } else if( inf_0 || inf_3 || - (orientation(p0, p1, p2, p3) == POSITIVE) ) { + } else if( inf_0 || inf_3 || (orientation(p0, p1, p2, p3) == POSITIVE) ) { // the two faces form a concavity, in which we might plug a tetrahedron int cospheric = 0; @@ -1825,14 +1816,14 @@ fill_hole_3D_ear( std::list & boundhole) bool on_unbounded_side = false; // we now look at all vertices that are on the boundary of the hole for(Surface::Vertex_iterator vit = surface.vertices_begin(); - (! violation ) && (vit != surface.vertices_end()); + (vit != surface.vertices_end()); vit++) { Vertex *v = (*vit).info(); if( (! is_infinite(Vertex_handle(v))) && (v != v0) && (v != v1) && (v != v2) && (v != v3)) { const Point & p = v->point(); - Bounded_side bs; // = side_of_sphere(v0, v1, v2, v3, p); + Bounded_side bs; if(inf_0) { bs = side_of_sphere_inf(p2, p1, p3, p); @@ -1841,195 +1832,163 @@ fill_hole_3D_ear( std::list & boundhole) } else { bs = Bounded_side( side_of_oriented_sphere(p0, p1, p2, p3, p) ); } - //CGAL_triangulation_assertion(bs == side_of_sphere(v0, v1, v2, v3, p)); + + if(bs == ON_BOUNDED_SIDE) { goto next_ear; } if((bs == ON_BOUNDARY)) { cospheric++; cospheric_vertices.insert(&(*vit)); } - violation = (bs == ON_BOUNDED_SIDE); - on_unbounded_side |= (bs == on_unbounded_side); + + on_unbounded_side |= (bs == ON_UNBOUNDED_SIDE); } } - + // we looked at all vertices // if there are cospheric points we have to test a little bit more - if( (! violation) && (cospheric > 0) ) { + if( cospheric > 0 ) { + std::set::iterator co_it, + not_found = cospheric_vertices.end(); if(inf_0 || inf_3) { - // the cospheric points are on the boundary of the convex hull we don't care + // the cospheric points are on the boundary of the convex hull if(! on_unbounded_side) { //std::cout << "all points are coplanar" << std::endl; } } else { - // for all edges that are incident to w2, check if the other vertex is cospheric - // and if the edge is coplanar to plane(v0, v2, v3) + // for all edges that are incident to w2, check if the other vertex + // is cospheric, and if the edge is coplanar to plane(v0, v2, v3) Vertex_circulator_3_2 vc = w2->incident_vertices(); Vertex_circulator_3_2 done = vc; - std::set::iterator co_it, not_found = cospheric_vertices.end(); + do { if( ! ((co_it = cospheric_vertices.find(&(*vc))) == not_found) ) { const Point & pc = (*co_it)->info()->point(); - violation = orientation(p0,p2,p3,pc) == COPLANAR; - } - } while( (! violation) && (++vc != done)); - if( ! violation ) { - // for all edges that are incident to w1, check if the other vertex is cospheric - // and if the edge is coplanar to plane(v0, v1, v3) - Vertex_circulator_3_2 vc = w1->incident_vertices(); - Vertex_circulator_3_2 done = vc; - std::set::iterator co_it, not_found = cospheric_vertices.end(); - do { - if( ! ((co_it = cospheric_vertices.find(&(*vc))) == not_found) ) { - const Point & pc = (*co_it)->info()->point(); - violation = (orientation(p0,p1,p3,pc) == COPLANAR); - } - } while( (! violation) && (++vc != done)); - } + if(orientation(p0,p2,p3,pc) == COPLANAR) { goto next_ear; } + } + } while( ++vc != done ); + + // for all edges that are incident to w1, check if the other vertex + // is cospheric, and if the edge is coplanar to plane(v0, v1, v3) + vc = w1->incident_vertices(); + done = vc; + do { + if( ! ((co_it = cospheric_vertices.find(&(*vc))) == not_found) ) { + const Point & pc = (*co_it)->info()->point(); + if(orientation(p0,p1,p3,pc) == COPLANAR) { goto next_ear; } + } + } while( ++vc != done ); + } - if( (! violation) && (! inf_0) && (! inf_1) && (! inf_2) && (! inf_3) ) { + if( (! inf_0) && (! inf_1) && (! inf_2) && (! inf_3) ) { + // all vertices are finite for(std::set::iterator it = cospheric_vertices.begin(); it != cospheric_vertices.end(); it++) { const Point & pit = (*it)->info()->point(); Vertex_circulator_3_2 vc = (*it)->incident_vertices(); Vertex_circulator_3_2 done = vc; - std::set::iterator co_it, not_found = cospheric_vertices.end(); do { - if( ! ((co_it = cospheric_vertices.find(&(*vc))) == not_found) ) { + if(! ((co_it = cospheric_vertices.find(&(*vc))) == not_found) ){ const Point & pc = (*co_it)->info()->point(); - violation = (orientation(p0,p3,pc, pit) == COPLANAR) - && (coplanar_orientation(p0,p3,pc, pit) == NEGATIVE); + if((orientation(p0,p3,pc, pit) == COPLANAR) + && (coplanar_orientation(p0,p3,pc, pit) == NEGATIVE)) { + goto next_ear; + } } - } while( (! violation) && (++vc != done)); + } while( ++vc != done); } } } // test a little bit more + // Face_3_2 *m_i = f->neighbor(i); Face_3_2 *m_j = f->neighbor(j); bool neighbor_i = m_i == n->neighbor(cw(fi)); bool neighbor_j = m_j == n->neighbor(ccw(fi)); - if((! violation) && (! (((! neighbor_i) && (! neighbor_j)) - && surface.is_edge(f->vertex(k), n->vertex(fi))))) { - // none of the vertices violates the Delaunay property - // and if we are in the flip case, the edge that would get introduced is not on the surface - // We are ready to plug the tetrahedron - // It may touch 2 triangles, + if( (((! neighbor_i) && (! neighbor_j)) + && surface.is_edge(f->vertex(k), n->vertex(fi)))) { + // We are in the flip case: + // The edge that would get introduced is on the surface + goto next_ear; + } + + // none of the vertices violates the Delaunay property + // We are ready to plug the tetrahedron - Cell_handle ch = create_cell(Vertex_handle(v0), Vertex_handle(v1), Vertex_handle(v2), Vertex_handle(v3), - NULL, NULL, NULL, NULL); - cells.push_back(ch); - Cell* c = handle2pointer(ch); - v0->set_cell(c); - v1->set_cell(c); - v2->set_cell(c); - v3->set_cell(c); - - //print(ch); - /* -Removal of point 265 : -7 -18 -15 + Cell_handle ch = create_cell(Vertex_handle(v0), + Vertex_handle(v1), + Vertex_handle(v2), + Vertex_handle(v3), + NULL, NULL, NULL, NULL); + cells.push_back(ch); + Cell* c = handle2pointer(ch); + v0->set_cell(c); + v1->set_cell(c); + v2->set_cell(c); + v3->set_cell(c); --9 -18 -14; -7 -21 -10; -8 -22 -16; -3 -21 -11; -flip --1 -20 -15; -8 -22 -16; -3 -21 -11; -9 -18 -14; -flip --9 -18 -14; -1 -20 -15; -8 -22 -16; -7 -19 -24; -remove --3 -11 -20; -7 -19 -24; -1 -20 -15; -9 -18 -14; -flip --9 -18 -14; -3 -11 -20; -7 -19 -24; -11 -14 -13; -remove --9 -18 -14; -3 -11 -20; -11 -14 -13; -7 -9 -11; -flip --11 -14 -13; -7 -9 -11; -5 -12 -9; -3 -11 -20; -flip --5 -12 -9; -9 -18 -14; -11 -14 -13; -7 -9 -11; -flip --5 -12 -9; -7 -21 -10; -9 -18 -14; -3 -21 -11; -remove --5 -12 -9; -3 -21 -11; -9 -18 -14; -1 -20 -15; -remove --3 -11 -20; -1 -20 -15; -5 -12 -9; -9 -18 -14; -remove --3 -11 -20; -9 -18 -14; -5 -12 -9; -7 -9 -11; -remove + // The new tetrahedron touches the faces that form the ear + Facet fac = n->info(); + c->set_neighbor(0, fac.first); + fac.first->set_neighbor(fac.second, c); + fac = f->info(); + c->set_neighbor(3, fac.first); + fac.first->set_neighbor(fac.second, c); -Unable to find an ear -4 4 2 --3 -11 -20; -5 -12 -9; -11 -14 -13; -7 -9 -11 - */ - - Facet fac = n->info(); - c->set_neighbor(0, fac.first); + // It may touch another face, + // or even two other faces if it is the last tetrahedron + if(neighbor_i) { + fac = m_i->info(); + c->set_neighbor(1, fac.first); fac.first->set_neighbor(fac.second, c); - fac = f->info(); - c->set_neighbor(3, fac.first); + } + if(neighbor_j) { + fac = m_j->info(); + c->set_neighbor(2, fac.first); fac.first->set_neighbor(fac.second, c); - CGAL_triangulation_assertion(c->neighbor(0) != c->neighbor(3)); - // 3, or even 4 if it is the last tetrahedron + } + + if((! neighbor_i) && (! neighbor_j)) { + surface.flip(f,k); + int ni = f->index(n); + f->set_edge(ni, false); + f->set_edge(cw(ni)); + f->set_edge(ccw(ni)); + last_op = f; k = -1; + + int fi = n->index(f); + n->set_edge(fi, false); + n->set_edge(cw(fi), f); + n->set_edge(ccw(fi), f); - if(neighbor_i) { - fac = m_i->info(); - c->set_neighbor(1, fac.first); - CGAL_triangulation_assertion(c->neighbor(1) != c->neighbor(3)); - CGAL_triangulation_assertion(c->neighbor(1) != c->neighbor(0)); - fac.first->set_neighbor(fac.second, c); - } - if(neighbor_j) { - fac = m_j->info(); - c->set_neighbor(2, fac.first); - CGAL_triangulation_assertion(c->neighbor(2) != c->neighbor(3)); - CGAL_triangulation_assertion(c->neighbor(2) != c->neighbor(0)); - if(neighbor_i) { - CGAL_triangulation_assertion(c->neighbor(2) != c->neighbor(1)); - } - fac.first->set_neighbor(fac.second, c); - } - - if((! neighbor_i) && (! neighbor_j)) { - surface.flip(f,k); - //std::cout << "flip" << std::endl; - int ni = f->index(n); - f->set_edge(ni, false); - f->set_edge(cw(ni)); - f->set_edge(ccw(ni)); - last_op = f; k = -1; - - int fi = n->index(f); - n->set_edge(fi, false); - n->set_edge(cw(fi), f); - n->set_edge(ccw(fi), f); - - f->set_info(std::make_pair(Cell_handle(c),2)); - n->set_info(std::make_pair(Cell_handle(c),1)); - } else if (neighbor_i && (! neighbor_j)) { - f->unlink(j); - surface.remove_degree_3(f->vertex(j), f); - //std::cout << "remove" << std::endl; - f->set_edge(); - last_op = f; k = -1; - f->set_info(std::make_pair(Cell_handle(c),2)); - } else if ((! neighbor_i) && neighbor_j) { - f->unlink(i); - surface.remove_degree_3(f->vertex(i), f); - //std::cout << "remove" << std::endl; - f->set_edge(); - last_op = f; k = -1; - f->set_info(std::make_pair(Cell_handle(c),1)); + f->set_info(std::make_pair(Cell_handle(c),2)); + n->set_info(std::make_pair(Cell_handle(c),1)); + } else if (neighbor_i && (! neighbor_j)) { + f->unlink(j); + surface.remove_degree_3(f->vertex(j), f); + f->set_edge(); + last_op = f; k = -1; + f->set_info(std::make_pair(Cell_handle(c),2)); + } else if ((! neighbor_i) && neighbor_j) { + f->unlink(i); + surface.remove_degree_3(f->vertex(i), f); + f->set_edge(); + last_op = f; k = -1; + f->set_info(std::make_pair(Cell_handle(c),1)); + } else { + if(surface.number_of_vertices() != 4) { + delete_cells(cells); + return false; } else { - if(surface.number_of_vertices() != 4) { - delete_cells(cells); - return false; - } else { - // when we leave the function the vertices and faces of the surface - // are deleted by the destructor - return true; - } + // when we leave the function the vertices and faces of the surface + // are deleted by the destructor + return true; } - }// if(! violation) - }// if(geom_traits().. + } + } // else if (inf_0 ... }// if(f->edge(k)) + next_ear: ; } // for(;;) }