From 2ae72afd4dfa76a5229e4e9685fe275cec20afa7 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 25 Feb 2010 15:29:19 +0000 Subject: [PATCH] Use a symbolic perturbation so that in the grid case all diagonals have the same orientation --- .../include/CGAL/Delaunay_triangulation_2.h | 10 +- .../include/CGAL/Regular_triangulation_2.h | 64 ++++++++++--- .../include/CGAL/Triangulation_2.h | 91 ++++++++++++++++++- .../CGAL/_test_cls_regular_triangulation_2.h | 2 +- .../include/CGAL/_test_cls_triangulation_2.h | 3 +- .../CGAL/_test_cls_triangulation_short_2.h | 7 +- 6 files changed, 151 insertions(+), 26 deletions(-) diff --git a/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h b/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h index f422dc73615..cf6d18d4f1f 100644 --- a/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h @@ -275,7 +275,7 @@ protected: int mi = this->_tds.mirror_index(f, i); Vertex_handle vm = this->_tds.mirror_vertex(f, i); if(this->is_infinite(vm)) continue; - if(this->side_of_oriented_circle(f, vm->point()) == ON_POSITIVE_SIDE) { + if(this->side_of_oriented_circle(f, vm->point(),true) == ON_POSITIVE_SIDE) { this->_tds.flip(f, i); edges.push_back(Edge(f, i)); edges.push_back(Edge(f, this->cw(i))); @@ -296,7 +296,7 @@ test_conflict(const Point &p, Face_handle fh) const // return true if P is inside the circumcircle of fh // if fh is infinite, return true when p is in the positive // halfspace or on the boundary and in the finite edge of fh - Oriented_side os = side_of_oriented_circle(fh,p); + Oriented_side os = side_of_oriented_circle(fh,p,true); if (os == ON_POSITIVE_SIDE) return true; if (os == ON_ORIENTED_BOUNDARY && is_infinite(fh)) { @@ -339,7 +339,7 @@ is_valid(bool verbose, int level) const for(int i=0; i<3; i++) { if ( ! is_infinite( this->mirror_vertex(it,i))) { result = result && ON_POSITIVE_SIDE != - side_of_oriented_circle( it, this->mirror_vertex(it,i)->point()); + side_of_oriented_circle( it, this->mirror_vertex(it,i)->point(), false); } CGAL_triangulation_assertion( result ); } @@ -420,7 +420,7 @@ look_nearest_neighbor(const Point& p, Vertex_handle& nn) const { Face_handle ni=f->neighbor(i); - if ( ON_POSITIVE_SIDE != side_of_oriented_circle(ni,p) ) return; + if ( ON_POSITIVE_SIDE != side_of_oriented_circle(ni,p,true) ) return; typename Geom_traits::Compare_distance_2 compare_distance = this->geom_traits().compare_distance_2_object(); @@ -576,7 +576,7 @@ propagating_flip(Face_handle& f,int i) Face_handle n = f->neighbor(i); if ( ON_POSITIVE_SIDE != - side_of_oriented_circle(n, f->vertex(i)->point()) ) { + side_of_oriented_circle(n, f->vertex(i)->point(), true) ) { return; } flip(f, i); diff --git a/Triangulation_2/include/CGAL/Regular_triangulation_2.h b/Triangulation_2/include/CGAL/Regular_triangulation_2.h index 8276ba8f61c..4af6fa60b4c 100644 --- a/Triangulation_2/include/CGAL/Regular_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Regular_triangulation_2.h @@ -25,6 +25,8 @@ #include #include +#include + CGAL_BEGIN_NAMESPACE template < typename K_ > @@ -192,14 +194,14 @@ public: Oriented_side power_test(const Weighted_point &p, const Weighted_point &q, const Weighted_point &r, - const Weighted_point &s) const; + const Weighted_point &s, bool perturb) const; Oriented_side power_test(const Weighted_point &p, const Weighted_point &q, const Weighted_point &r) const; Oriented_side power_test(const Weighted_point &p, const Weighted_point &r) const; Oriented_side power_test(const Face_handle &f, - const Weighted_point &p) const; + const Weighted_point &p, bool perturb=false) const; Oriented_side power_test(const Face_handle& f, int i, const Weighted_point &p) const; @@ -615,7 +617,7 @@ operator=(const Self &tr) template < class Gt, class Tds > Oriented_side Regular_triangulation_2:: -power_test(const Face_handle &f, const Weighted_point &p) const +power_test(const Face_handle &f, const Weighted_point &p, bool perturb) const { // p is supposed to be a finite point // if f is a finite face, @@ -627,7 +629,7 @@ power_test(const Face_handle &f, const Weighted_point &p) const if ( ! f->has_vertex(infinite_vertex(), i) ) return power_test(f->vertex(0)->point(), f->vertex(1)->point(), - f->vertex(2)->point(),p); + f->vertex(2)->point(),p, perturb); Orientation o = orientation(f->vertex(ccw(i))->point(), f->vertex( cw(i))->point(), @@ -662,12 +664,48 @@ template < class Gt, class Tds > inline Oriented_side Regular_triangulation_2:: -power_test(const Weighted_point &p, - const Weighted_point &q, - const Weighted_point &r, - const Weighted_point &s) const +power_test(const Weighted_point &p0, + const Weighted_point &p1, + const Weighted_point &p2, + const Weighted_point &p, + bool perturb) const { - return geom_traits().power_test_2_object()(p,q,r,s); + CGAL_triangulation_precondition( orientation(p0, p1, p2) == POSITIVE ); + + using namespace boost; + + Oriented_side os = geom_traits().power_test_2_object()(p0, p1, p2, p); + + if ( (os != ON_ORIENTED_BOUNDARY) || (! perturb)) + return os; + + // We are now in a degenerate case => we do a symbolic perturbation. + + // We sort the points lexicographically. + const Weighted_point * points[4] = {&p0, &p1, &p2, &p}; + std::sort(points, points + 4, + bind(geom_traits().compare_xy_2_object(), + bind(Dereference(), _1), + bind(Dereference(), _2)) == SMALLER); + + // We successively look whether the leading monomial, then 2nd monomial + // of the determinant has non null coefficient. + // 2 iterations are enough (cf paper) + for (int i=3; i>1; --i) { + if (points[i] == &p) + return ON_NEGATIVE_SIDE; // since p0 p1 p2 are non collinear + // and positively oriented + Orientation o; + if (points[i] == &p2 && (o = orientation(p0,p1,p)) != COLLINEAR ) + return o; + if (points[i] == &p1 && (o = orientation(p0,p,p2)) != COLLINEAR ) + return o; + if (points[i] == &p0 && (o = orientation(p,p1,p2)) != COLLINEAR ) + return o; + } + + CGAL_triangulation_assertion(false); + return ON_NEGATIVE_SIDE; } template < class Gt, class Tds > @@ -1090,7 +1128,7 @@ insert(const Weighted_point &p, Locate_type lt, Face_handle loc, int li) { CGAL_precondition (dimension() >= 1); Oriented_side os = dimension() == 1 ? power_test (loc, li, p) : - power_test (loc, p); + power_test (loc, p, true); if (os < 0) { if (is_infinite (loc)) loc = loc->neighbor (li); @@ -1103,7 +1141,7 @@ insert(const Weighted_point &p, Locate_type lt, Face_handle loc, int li) case Base::FACE: { CGAL_precondition (dimension() >= 2); - if (power_test (loc, p) < 0) { + if (power_test (loc, p, true) < 0) { return hide_new_vertex(loc,p); } v = insert_in_face (p, loc); @@ -1514,7 +1552,7 @@ fill_hole_regular(std::list & first_hole) p2=p; cut_after=hit; } - else if (power_test(p0,p1,p2,p) == + else if (power_test(p0,p1,p2,p,true) == ON_POSITIVE_SIDE) { v2=vv; @@ -1824,7 +1862,7 @@ stack_flip(Vertex_handle v, Faces_around_stack &faces_around) // now dimension() == 2 //test the regularity of edge (f,i) //if( power_test(n, v->point()) == ON_NEGATIVE_SIDE) - if( power_test(n, v->point()) != ON_POSITIVE_SIDE) + if( power_test(n, v->point(), true) != ON_POSITIVE_SIDE) return; if(is_infinite(f,i)) diff --git a/Triangulation_2/include/CGAL/Triangulation_2.h b/Triangulation_2/include/CGAL/Triangulation_2.h index 37891ea6a9c..0e0ed813c3d 100644 --- a/Triangulation_2/include/CGAL/Triangulation_2.h +++ b/Triangulation_2/include/CGAL/Triangulation_2.h @@ -91,6 +91,20 @@ public: typedef typename Tds::Vertex_iterator All_vertices_iterator; + class Perturbation_order { + const Self *t; + + public: + Perturbation_order(const Self *tr) + : t(tr) {} + + bool operator()(const Point *p, const Point *q) const { + return t->compare_xy(*p, *q) == SMALLER; + } + }; + + friend class Perturbation_order; + // This class is used to generate the Finite_*_iterators. class Infinite_tester { @@ -363,14 +377,21 @@ public: Oriented_side oriented_side(Face_handle f, const Point &p) const; + +Oriented_side +side_of_oriented_circle(const Point &p0, const Point &p1, const Point &p2, + const Point &p, bool perturb) const; + Oriented_side - side_of_oriented_circle(Face_handle f, const Point & p) const; + side_of_oriented_circle(Face_handle f, const Point & p, bool perturb = false) const; bool collinear_between(const Point& p, const Point& q, const Point& r) const; Comparison_result compare_x(const Point& p, const Point& q) const; + + Comparison_result compare_xy(const Point& p, const Point& q) const; Comparison_result compare_y(const Point& p, const Point& q) const; bool xy_equal(const Point& p, const Point& q) const; Orientation orientation(const Point& p, @@ -1459,7 +1480,8 @@ fill_hole_delaunay(std::list & first_hole) if (orientation(p0,p1,p) == COUNTERCLOCKWISE) { if (is_infinite(v2)) { v2=vv; v3=vv; cut_after=hit;} else{ - if (in_circle(p0,p1,v3->point(),p) == ON_POSITIVE_SIDE){ + // + if (this->side_of_oriented_circle(p0,p1,v3->point(),p,true) == ON_POSITIVE_SIDE){ v2=vv; v3=vv; cut_after=hit;} } } @@ -2432,17 +2454,67 @@ oriented_side(Face_handle f, const Point &p) const p); } + +template +Oriented_side +Triangulation_2:: +side_of_oriented_circle(const Point &p0, const Point &p1, const Point &p2, + const Point &p, bool perturb) const +{ + CGAL_triangulation_precondition( orientation(p0, p1, p2) == POSITIVE ); + + Gt::Side_of_oriented_circle_2 pred = geom_traits().side_of_oriented_circle_2_object(); + Oriented_side os = + pred(p0, p1, p2, p); + if ((os != ON_ORIENTED_BOUNDARY) || (! perturb)) + return os; + + // We are now in a degenerate case => we do a symbolic perturbation. + + // We sort the points lexicographically. + const Point * points[4] = {&p0, &p1, &p2, &p}; + std::sort(points, points+4, Perturbation_order(this) ); + + // We successively look whether the leading monomial, then 2nd monomial + // of the determinant has non null coefficient. + // 2 iterations are enough (cf paper) + for (int i=3; i>0; --i) { + if (points[i] == &p) + return ON_NEGATIVE_SIDE; // since p0 p1 p2 are non collinear + // and positively oriented + Orientation o; + if (points[i] == &p2 && (o = orientation(p0,p1,p)) != COLLINEAR ) + return Oriented_side(o); + if (points[i] == &p1 && (o = orientation(p0,p,p2)) != COLLINEAR ) + return Oriented_side(o); + if (points[i] == &p0 && (o = orientation(p,p1,p2)) != COLLINEAR ) + return Oriented_side(o); + } + CGAL_triangulation_assertion(false); + return ON_NEGATIVE_SIDE; +} + + + + + + template < class Gt, class Tds > Oriented_side Triangulation_2:: -side_of_oriented_circle(Face_handle f, const Point & p) const +side_of_oriented_circle(Face_handle f, const Point & p, bool perturb) const { if ( ! is_infinite(f) ) { + /* typename Gt::Side_of_oriented_circle_2 in_circle = geom_traits().side_of_oriented_circle_2_object(); return in_circle(f->vertex(0)->point(), f->vertex(1)->point(), f->vertex(2)->point(),p); + */ + return this->side_of_oriented_circle(f->vertex(0)->point(), + f->vertex(1)->point(), + f->vertex(2)->point(),p, perturb); } int i = f->index(infinite_vertex()); @@ -2488,6 +2560,19 @@ compare_x(const Point& p, const Point& q) const return geom_traits().compare_x_2_object()(p,q); } +template +inline +Comparison_result +Triangulation_2:: +compare_xy(const Point& p, const Point& q) const +{ + Comparison_result res = geom_traits().compare_x_2_object()(p,q); + if(res == EQUAL){ + return geom_traits().compare_y_2_object()(p,q); + } + return res; +} + template inline Comparison_result diff --git a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_regular_triangulation_2.h b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_regular_triangulation_2.h index 57899088a9b..30d4437b4e5 100644 --- a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_regular_triangulation_2.h +++ b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_regular_triangulation_2.h @@ -555,7 +555,7 @@ _test_cls_regular_triangulation_2( const Triangulation & ) assert(fc==fc2); n=0; do {fc2++ ; n = n+1;} while (fc2 != fc); - assert(n==3); + assert(n==4); /*****************************/ /******** Miscellaneaous *****/ diff --git a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_2.h b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_2.h index b2393233be1..b5b44678b84 100644 --- a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_2.h +++ b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_2.h @@ -311,7 +311,8 @@ _test_cls_triangulation_2( const Triangul & ) assert(!T2_8.is_infinite(ff)); Face_handle f2 = ff->neighbor(li); assert(!T2_8.is_infinite(f2)); - T2_8.flip(ff,0); + int fli = ff->index(f2); + T2_8.flip(ff,fli); assert( T2_8.is_valid() ); //make_hole star_hole diff --git a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_short_2.h b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_short_2.h index 05be9f7729d..60ddf052d35 100644 --- a/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_short_2.h +++ b/Triangulation_2/test/Triangulation_2/include/CGAL/_test_cls_triangulation_short_2.h @@ -234,7 +234,8 @@ _test_cls_triangulation_short_2( const Triangul &) assert(!T2_8.is_infinite(ff)); Face_handle f2 = ff->neighbor(li); assert(!T2_8.is_infinite(f2)); - T2_8.flip(ff,0); + int fli = ff->index(f2); + T2_8.flip(ff,fli); assert( T2_8.is_valid() ); @@ -412,7 +413,7 @@ _test_cls_triangulation_short_2( const Triangul &) do {fc2++ ; n = n+1;} while (fc2 != fc); assert(T2_8.number_of_vertices()>=2); assert(T2_8.is_valid()); - fc= T2_8.line_walk(Point(5,4,10),Point(5,5)); + fc= T2_8.line_walk(Point(1,4,10),Point(20,5,10)); fc2=fc; n=0; assert(fc==fc2); @@ -431,7 +432,7 @@ _test_cls_triangulation_short_2( const Triangul &) assert(fc==fc2); n=0; do {fc2++ ; n = n+1;} while (fc2 != fc); - assert(n==3); + assert(n==4 || n==3); /*****************************/ /******** Miscellaneaous *****/