From c1959465eb3c56c7bdee3ebca37cb39b5c2de2bf Mon Sep 17 00:00:00 2001 From: Monique Teillaud Date: Mon, 19 Jul 1999 14:59:52 +0000 Subject: [PATCH] INITIAL VERSION ! only for a triangulation where the first four points give a triangulation with dimension 3 seems to work... --- .../include/CGAL/Regular_triangulation_3.h | 668 +++++++----------- ...Regular_triangulation_euclidean_traits_3.h | 7 +- 2 files changed, 254 insertions(+), 421 deletions(-) diff --git a/Packages/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Packages/Triangulation_3/include/CGAL/Regular_triangulation_3.h index c600b7ac484..b29b7c259a5 100644 --- a/Packages/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Packages/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -14,9 +14,10 @@ // file : include/CGAL/Regular_triangulation_3.h // revision : $Revision$ // revision_date : $Date$ -// author(s) : Sophie Balaven +// author(s) : Monique Teillaud +// Sophie Balaven // Sylvain Pion -// Monique Teillaud +// // // coordinator : Mariette Yvinec // @@ -29,10 +30,6 @@ #include #include -#define DEBUG_INSERT 0 -#define DEBUG_STAR 0 -#define DEBUG_SIDE_OF_SPHERE 0 - CGAL_BEGIN_NAMESPACE template < class Gt, class Tds> @@ -120,434 +117,269 @@ public: } #endif // CGAL_CFG_NO_MEMBER_TEMPLATES - /// +++ Insertion +++ // + Vertex_handle insert( const Weighted_point &p ); - // -- 1 -- // - Vertex_handle insert( const Weighted_point &p ) - { - Cell_handle start; - if ( dimension() >= 1 ) { - // there is at least one finite "cell" (or facet or edge) - start = infinite_vertex()->cell() - ->neighbor(infinite_vertex()->cell()->index(infinite_vertex())); - } - else { - start = NULL; - } - return insert( p, start ); - } - - // -- 2 -- // - Vertex_handle insert( const Weighted_point & p, Cell_handle start ) { - - Vertex_handle v; + Vertex_handle insert( const Weighted_point & p, Cell_handle start ); - if (DEBUG_INSERT){ - std::cout << "dimension dela triangulation avant d'inserer le nouveau point : "; - } - - switch (dimension()) { - case 3: - { - if (DEBUG_INSERT){ - std::cout << "3" << std::endl; - } - - Locate_type lt; - int li, lj; - Cell_handle c = locate( p, start, lt, li, lj); - - switch (lt) { - case OUTSIDE_CONVEX_HULL: - case CELL: - case FACET: - case EDGE: - { - if (DEBUG_INSERT){ - std::cout << "j'essaye d'inserer un nouveau point en dimension 3\n"; - std::cout << "-- dim = 3 -- " << p << std::endl; - } - v = new Vertex(p); - if (DEBUG_INSERT){ - std::cout << "vertex : " << v->point() << std::endl; - } - // a voir : comment compter les sommets redondants ? (***) - // set_number_of_vertices(number_of_vertices()+1); - std::set > conflicts; - std::set > deleted_points; - Cell_handle aconflict; - int ineighbor; - if (in_conflict_3(p, c)) { - find_conflicts_3(conflicts, p, c, aconflict, ineighbor); - deleted_points = - star_region_delete_points(conflicts,&(*v),&(*aconflict), - ineighbor); - // debug - if (DEBUG_INSERT){ - std::set >::const_iterator it; - std::cout << "##### Points supprimes #####\n"; - for( it = deleted_points.begin(); it != deleted_points.end(); ++it) { - Point *p_tmp = (Point *) *it; - std::cout << p_tmp->point() << std::endl; - } - std::cout << "############################\n"; - } - - // a voir : comment compter les sommets redondants ? (***) - set_number_of_vertices(number_of_vertices()+1); - } - // else : traiter le cas des points redondants a stocker dans - // la face associee pour le cas d'une future suppression - return v; - } - case VERTEX: - return c->vertex(li); - default : - CGAL_triangulation_assertion(false); // impossible - } - break; - } - case 2: - { - if (DEBUG_INSERT){ - std::cout << "2" << std::endl; - } - Locate_type lt; - int li, lj; - Cell_handle c = locate( p, start, lt, li, lj); - switch (lt) { - case OUTSIDE_CONVEX_HULL: - case CELL: - case FACET: - case EDGE: - { - // 2D : particular case - // not yet implemented - } - case VERTEX: - return c->vertex(li); - case OUTSIDE_AFFINE_HULL: - { - // if the 2d triangulation is Regular, the 3d - // triangulation will be Regular - if (DEBUG_INSERT){ - std::cout << "j'insere un point en dehors de l'enveloppe affine\n"; - std::cout << "-- dim = 2 -- " << p << std::endl; - } - return - Triangulation_3::insert_outside_affine_hull(p); - } - } - } - default : // dim <= 1 - if (DEBUG_INSERT){ - if (dimension() == 1) std::cout << "1" << std::endl; - else std::cout << "0" << std::endl; - std::cout << "j'insere un point dans la triangulation de base\n"; - std::cout << "-- dim <= 1 -- " << p << std::endl; - } - return Triangulation_3::insert(p); - } - return Triangulation_3::insert(p); // to avoid warning with egcs - } + Bounded_side + side_of_power_sphere( Cell_handle c, const Weighted_point &p) const; + bool is_valid(bool verbose = false, int level = 0) const; + +private: + bool in_conflict_3(const Weighted_point & p, const Cell_handle & c) + { + return side_of_power_sphere(c, p) == ON_BOUNDED_SIDE; + } + + void find_conflicts_3(std::set > &conflicts, + const Weighted_point & p, + Cell_handle c, Cell_handle & ac, int & i); std::set > - star_region_delete_points( std::set > & region, Vertex* v, - Cell* c, int li) + star_region_delete_points( std::set > & region, + Vertex* v, + Cell* c, int li); // region is a set of connected cells // c belongs to region and has facet i on the boundary of region // replaces the cells in region // by linking v to the boundary of region // deleted weighted points that are not in the triangulation // anymore, pts will be the list of deleted points - { - std::set > vert; - Cell *c_tmp; - Vertex *v_tmp; - Vertex_handle vh; - std::set > pts; - Point *p; - int i=0; - - if (DEBUG_STAR) - std::cout << "@ star_region_delete_points\n"; - // for each cell to be deleted, keep vertices - std::set >::const_iterator it; - for( it = region.begin(); it != region.end(); ++it) { - c_tmp = (Cell *) *it; - if (DEBUG_STAR) { - std::cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"; - std::cout << "Region a detruire : cellule " << i++ << std::endl; - } - // store vertices - for (int i=0; i<4 ; i++){ - vh = c_tmp->vertex(i); - if (DEBUG_STAR) - std::cout << *vh << std::endl; - if ( (vert.find((void*) &(*(vh))) == vert.end()) && - (! is_infinite(v_tmp)) ) { - if (DEBUG_STAR) - std::cout << "stocke le sommet : " << *vh << std::endl; - vert.insert( (void*) &( *(vh) ) ); - } - } - } - if (DEBUG_STAR) - std::cout << "vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv\n"; - - - if (DEBUG_STAR) - std::cout << "j'appelle star_region\n"; - - // Create the new faces and delete old ones - _tds.star_region( region, v, c, li ); - - // get the vertices incident to v - std::set > inc_vert; - incident_vertices(v, inc_vert); - - // for each vertex, check if it is a vertex incident to v - // if not, delete it - if (DEBUG_STAR) - std::cout << "je passe en revue les sommets\n"; - for( it = vert.begin(); it != vert.end(); ++it) { - v_tmp = (Vertex *) *it; - if (DEBUG_STAR) - std::cout << *v_tmp << std::endl; - if ( (inc_vert.find( v_tmp )) == inc_vert.end() ) { - // vertex has to be deleted and point to be stored - p = new Point( v_tmp->point() ); - if (DEBUG_STAR) { - std::cout << "-- Suppression du point -- " << *p << std::endl; - } - pts.insert( p ); - set_number_of_vertices(number_of_vertices()-1); - delete(*it); - } - } - - // returns list of deleted points - return pts; - } - - - -private: - bool in_conflict_3(const Weighted_point & p, const Cell_handle & c) - { - return side_of_sphere(c, p) == ON_BOUNDED_SIDE; - } - - void find_conflicts_3(std::set > &conflicts, - const Weighted_point & p, - Cell_handle c, Cell_handle & ac, int & i) - { - if ( ( conflicts.find( (void *) &(*c) ) ) != conflicts.end() ) - { - return; // c was already found - } - - (void) conflicts.insert( (void *) &(*c) ); - - for ( int j=0; j<4; j++ ) { - if ( in_conflict_3( p, c->neighbor(j) ) ) { - find_conflicts_3(conflicts, p, c->neighbor(j), ac, i); - } - else { - ac = c; - i = j; - } - } - } - - - void find_conflicts_2(std::set > &conflicts, - const Weighted_point & p, - Cell_handle c, Cell_handle & ac, int & i) - { - // find_conflicts_2 : not yet implemented - ; - } - -public: - - Bounded_side - side_of_sphere( Cell_handle c, const Weighted_point &p) const - { - CGAL_triangulation_precondition( dimension() == 3 ); - int i3; - if ( ! c->has_vertex( infinite_vertex(), i3 ) ) { - // Cas d'une face de dimension finie - Oriented_side o = geom_traits().power_test_3(c->vertex(0)->point(), - c->vertex(1)->point(), - c->vertex(2)->point(), - c->vertex(3)->point(),p); - return Bounded_side(o); - } - else { - // cas d'une face infinie - int i0,i1,i2; - if ( (i3%2) == 1 ) { - i0 = (i3+1)&3; - i1 = (i3+2)&3; - i2 = (i3+3)&3; - } - else { - i0 = (i3+2)&3; - i1 = (i3+1)&3; - i2 = (i3+3)&3; - } - - bool coll = false; - Oriented_side side; - // on verifie que p n'est pas collineaire a 2 autres points - // dans ce cas, on traite un cas en dimension 1 - if (geom_traits().collinear(c->vertex(i0)->point(), - c->vertex(i1)->point(), p)) - { - if (DEBUG_SIDE_OF_SPHERE) { - std::cout << "le point est collineaire a i0("; - std::cout << c->vertex(i0)->point() << ") et i1("; - std::cout << c->vertex(i1)->point() << ")\n"; - } - coll = true; - side = geom_traits().power_test_3 - ( c->vertex(i0)->point(), c->vertex(i1)->point(), p ); - } - else if (geom_traits().collinear(c->vertex(i0)->point(), - c->vertex(i2)->point(), p)) - { - if (DEBUG_SIDE_OF_SPHERE) { - std::cout << "le point est collineaire a i0("; - std::cout << c->vertex(i0)->point() << ") et i2("; - std::cout << c->vertex(i2)->point() << ")\n"; - } - coll = true; - side = geom_traits().power_test_3 - ( c->vertex(i0)->point(), c->vertex(i2)->point(), p ); - } - else if (geom_traits().collinear(c->vertex(i1)->point(), - c->vertex(i2)->point(), p)) - { - if (DEBUG_SIDE_OF_SPHERE) { - std::cout << "le point est collineaire a i1("; - std::cout << c->vertex(i1)->point() << ") et i2("; - std::cout << c->vertex(i2)->point() << ")\n"; - } - coll = true; - side = geom_traits().power_test_3 - ( c->vertex(i1)->point(), c->vertex(i2)->point(), p ); - } - - if (coll) - return Bounded_side(o); - - // teste de quel cote du plan defini par (i0, i1, i2) p repose - Orientation - o = geom_traits().orientation(c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), - p); - - if (o != ZERO) - return Bounded_side(o); - - // i0, i1, i2, p - coplanar - //std::cout << "les points sont coplanaires\n"; - Oriented_side s = - geom_traits().power_test_3 - ( c->vertex(i0)->point(), - c->vertex(i1)->point(), - c->vertex(i2)->point(), - p ); - - return Bounded_side(s); - } - return ON_UNBOUNDED_SIDE;// to avoid warning with egcs - }// end side of sphere - - - - - Bounded_side side_of_circle( const Facet &f, const Weighted_point &p) const; - Bounded_side side_of_circle( Cell_handle c, int i, - const Weighted_point &p) const; - - bool is_valid(bool verbose = false, int level = 0) const - { - if ( ! tds().is_valid(verbose,level) ) { - if (verbose) { std::cerr << "invalid data structure" << std::endl; } - CGAL_triangulation_assertion(false); return false; - } - if ( &(*infinite_vertex()) == NULL ) { - if (verbose) { std::cerr << "no infinite vertex" << std::endl; } - CGAL_triangulation_assertion(false); return false; - } - - int i; - - switch ( dimension() ) { - case 3: - { - Cell_iterator it; - for ( it = finite_cells_begin(); it != cells_end(); ++it ) { - is_valid_finite((*it).handle()); - for ( i=0; i<4; i++ ) { - if ( side_of_sphere - ( (*it).handle(), - it->vertex( (it->neighbor(i))->index((*it).handle() ) ) - ->point() ) - == ON_BOUNDED_SIDE ) { - if (verbose) { - std::cerr << "non-empty sphere " << std::endl; - } - CGAL_triangulation_assertion(false); return false; - } - } - } - break; - } - // case 2: - // { - // Facet_iterator it; - // for ( it = finite_facets_begin(); it != facets_end(); ++it ) { - // is_valid_finite((*it).first); - // for ( i=0; i<2; i++ ) { - // if ( side_of_circle - // ( (*it).first, 3, - // (*it).first - // ->vertex( (((*it).first)->neighbor(i))->index((*it).first) )->point() ) - // == ON_BOUNDED_SIDE ) { - // if (verbose) { - // std::cerr << "non-empty circle " << std::endl; - // } - // CGAL_triangulation_assertion(false); return false; - // } - // } - // } - // break; - // } - case 1: - { - Edge_iterator it; - for ( it = finite_edges_begin(); it != edges_end(); ++it ) { - is_valid_finite((*it).first); - } - break; - } - } - if (verbose) { std::cerr << "Regular valid triangulation" << std::endl;} - return true; - } }; -/////////////////////////// IMPLEMENTATION ////////////////////////// +template < class Gt, class Tds > +void +Regular_triangulation_3:: +find_conflicts_3(std::set > &conflicts, + const Weighted_point & p, + Cell_handle c, Cell_handle & ac, int & i) + // DUPLICATED CODE !!! see Delaunay +{ + if ( ( conflicts.find( (void *) &(*c) ) ) != conflicts.end() ) + return; // c was already found + + (void) conflicts.insert( (void *) &(*c) ); + + for ( int j=0; j<4; j++ ) { + if ( in_conflict_3( p, c->neighbor(j) ) ) { + find_conflicts_3(conflicts, p, c->neighbor(j), ac, i); + } + else { + ac = c; + i = j; + } + } +} -// +++ Insertion +++ // +template < class Gt, class Tds > +Bounded_side +Regular_triangulation_3:: +side_of_power_sphere( Cell_handle c, const Weighted_point &p) const +{ + CGAL_triangulation_precondition( dimension() == 3 ); + int i3; + if ( ! c->has_vertex( infinite_vertex(), i3 ) ) { + return Bounded_side( geom_traits().power_test + (c->vertex(0)->point(), + c->vertex(1)->point(), + c->vertex(2)->point(), + c->vertex(3)->point(),p) ); + } + // else infinite cell : + int i0,i1,i2; + if ( (i3%2) == 1 ) { + i0 = (i3+1)&3; + i1 = (i3+2)&3; + i2 = (i3+3)&3; + } + else { + i0 = (i3+2)&3; + i1 = (i3+1)&3; + i2 = (i3+3)&3; + } -// +++ Determination de la region en conflit +++ // + // general case + Orientation + o = geom_traits().orientation(c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), + p); + if (o != ZERO) + return Bounded_side(o); + // else p coplanar with i0,i1,i2 + return Bounded_side( geom_traits().power_test + ( c->vertex(i0)->point(), + c->vertex(i1)->point(), + c->vertex(i2)->point(), p ) ); +}// end side of power sphere + +template < class Gt, class Tds > +Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert(const Weighted_point & p ) +{ + Cell_handle start; + if ( dimension() >= 1 ) { + // there is at least one finite "cell" (or facet or edge) + start = infinite_vertex()->cell() + ->neighbor(infinite_vertex()->cell()->index(infinite_vertex())); + } + else { + start = NULL; + } + return Regular_triangulation_3::insert( p, start ); +} + +template < class Gt, class Tds > +Regular_triangulation_3::Vertex_handle +Regular_triangulation_3:: +insert(const Weighted_point & p, Cell_handle start ) +{ + Vertex_handle v; + + switch (dimension()) { + case 3: + { + Locate_type lt; + int li, lj; + Cell_handle c = locate( p, start, lt, li, lj); + + if ( lt == VERTEX ) { return c->vertex(li); } + // TBD : look at the weight... + + // else + v = new Vertex(p); + std::set > conflicts; + std::set > deleted_points; + Cell_handle aconflict; + int ineighbor; + if (in_conflict_3(p, c)) { + find_conflicts_3(conflicts, p, c, aconflict, ineighbor); + deleted_points = + star_region_delete_points(conflicts,&(*v),&(*aconflict), + ineighbor); + // a voir : comment compter les sommets redondants ? (***) + set_number_of_vertices(number_of_vertices()+1); + } + // else : traiter le cas des points redondants a stocker dans + // la face associee pour le cas d'une future suppression + return v; + } + default : + { + // temporary : will only work for non degenerated dimensions + // (only for the first 4 points) + return Triangulation_3::insert(p,start); + } + } +} + +template < class Gt, class Tds > +std::set > +Regular_triangulation_3:: +star_region_delete_points(std::set > & region, + Vertex* v, Cell* c, int li) + // region is a set of connected cells + // c belongs to region and has facet i on the boundary of region + // replaces the cells in region + // by linking v to the boundary of region + // deleted weighted points that are not in the triangulation + // anymore, pts will be the list of deleted points +{ + std::set > vert; + Cell *c_tmp; + Vertex *v_tmp; + Vertex_handle vh; + std::set > pts; + Point *p; + + // for each cell to be deleted, keep vertices + std::set >::const_iterator it; + for( it = region.begin(); it != region.end(); ++it) { + c_tmp = (Cell *) *it; + // store vertices + for (int i=0; i<4 ; i++){ + vh = c_tmp->vertex(i); + if ( (vert.find((void*) &(*(vh))) == vert.end()) && + (! is_infinite(vh)) ) { + vert.insert( (void*) &(*(vh) ) ); + } + } + } + + // Create the new faces and delete old ones + _tds.star_region( region, v, c, li ); + + // get the vertices incident to v + std::set > inc_vert; + incident_vertices(v, inc_vert); + + // for each vertex, check if it is a vertex incident to v + // if not, delete it + for( it = vert.begin(); it != vert.end(); ++it) { + v_tmp = (Vertex *) *it; + if ( (inc_vert.find( v_tmp )) == inc_vert.end() ) { + // vertex has to be deleted and point to be stored + p = new Point( v_tmp->point() ); + pts.insert( p ); + set_number_of_vertices(number_of_vertices()-1); + delete(*it); + } + } + + // returns list of deleted points + return pts; +} + +template < class Gt, class Tds > +bool +Regular_triangulation_3:: +is_valid(bool verbose = false, int level = 0) const +{ + if ( ! tds().is_valid(verbose,level) ) { + if (verbose) { std::cerr << "invalid data structure" << std::endl; } + CGAL_triangulation_assertion(false); return false; + } + + if ( &(*infinite_vertex()) == NULL ) { + if (verbose) { std::cerr << "no infinite vertex" << std::endl; } + CGAL_triangulation_assertion(false); return false; + } + + int i; + switch ( dimension() ) { + case 3: + { + Cell_iterator it; + for ( it = finite_cells_begin(); it != cells_end(); ++it ) { + is_valid_finite((*it).handle()); + for ( i=0; i<4; i++ ) { + if ( side_of_power_sphere + ( (*it).handle(), + it->vertex( (it->neighbor(i))->index((*it).handle() ) ) + ->point() ) + == ON_BOUNDED_SIDE ) { + if (verbose) { + std::cerr << "non-empty sphere " << std::endl; + } + CGAL_triangulation_assertion(false); return false; + } + } + } + break; + } + default: + { + std::cerr << "degenerate dimensions not yet implemented" + << std::endl; + CGAL_triangulation_assertion(false); + } + } + if (verbose) { std::cerr << "Regular valid triangulation" << std::endl;} + return true; +} CGAL_END_NAMESPACE #endif CGAL_REGULAR_TRIANGULATION_3_H diff --git a/Packages/Triangulation_3/include/CGAL/Regular_triangulation_euclidean_traits_3.h b/Packages/Triangulation_3/include/CGAL/Regular_triangulation_euclidean_traits_3.h index a4fd52b4a53..bfece7d7c1b 100644 --- a/Packages/Triangulation_3/include/CGAL/Regular_triangulation_euclidean_traits_3.h +++ b/Packages/Triangulation_3/include/CGAL/Regular_triangulation_euclidean_traits_3.h @@ -42,14 +42,14 @@ CGAL_BEGIN_NAMESPACE -template < class R, class W = typename R::FT> +template < class R, class W = typename R::RT> class Regular_triangulation_euclidean_traits_3 - : public Triangulation_euclidean_traits_3 + : public Triangulation_geom_traits_3 { public: typedef R Rep; typedef W Weight; - typedef Triangulation_euclidean_traits_3 Traits; + typedef Triangulation_geom_traits_3 Traits; typedef Traits::Point Bare_point; typedef Weighted_point Weighted_point; typedef Weighted_point Point; @@ -72,6 +72,7 @@ public: const Weighted_point &t) const { CGAL_triangulation_precondition( ! collinear(p, q, r) ); + CGAL_triangulation_precondition( coplanar(p,q,r,t) ); return CGAL::power_test(p, q, r, t); }