Modifications to code (faster insertion via spatial sorting) + new test for Euclidean DT_2

This commit is contained in:
Iordan Iordanov 2017-02-28 15:51:39 +01:00
parent fef53805d5
commit cab6689327
6 changed files with 407 additions and 318 deletions

View File

@ -49,158 +49,167 @@ namespace CGAL {
>
class Periodic_4_hyperbolic_Delaunay_triangulation_2: public Periodic_4_hyperbolic_triangulation_2<GT, TDS> {
typedef Periodic_4_hyperbolic_Delaunay_triangulation_2<GT, TDS> Self;
typedef Periodic_4_hyperbolic_triangulation_2<GT, TDS> Base;
typedef Periodic_4_hyperbolic_Delaunay_triangulation_2<GT, TDS> Self;
typedef Periodic_4_hyperbolic_triangulation_2<GT, TDS> Base;
public:
#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2
using Base::cw;
using Base::ccw;
using Base::geom_traits;
using Base::tds;
using Base::cw;
using Base::ccw;
using Base::geom_traits;
using Base::tds;
#endif
typedef typename Base::Locate_type Locate_type;
typedef typename Base::Geometric_traits Geometric_traits;
typedef typename Base::Triangulation_data_structure Triangulation_data_structure;
//typedef typename Base::Int Int;
typedef typename Base::Offset Offset;
typedef typename Base::Circle_2 Circle_2;
typedef Circle_2 Circle;
typedef typename Base::Point_2 Point_2;
typedef Point_2 Point;
typedef typename Base::Segment_2 Segment_2;
typedef Segment_2 Segment;
typedef typename Base::Triangle_2 Triangle_2;
typedef Triangle_2 Triangle;
typedef typename Base::Locate_type Locate_type;
typedef typename Base::Geometric_traits Geometric_traits;
typedef typename Base::Triangulation_data_structure Triangulation_data_structure;
//typedef typename Base::Int Int;
typedef typename Base::Offset Offset;
typedef typename Base::Circle_2 Circle_2;
typedef Circle_2 Circle;
typedef typename Base::Point_2 Point_2;
typedef Point_2 Point;
typedef typename Base::Segment_2 Segment_2;
typedef Segment_2 Segment;
typedef typename Base::Triangle_2 Triangle_2;
typedef Triangle_2 Triangle;
typedef typename Base::Periodic_point Periodic_point;
typedef typename Base::Periodic_segment Periodic_segment;
typedef typename Base::Periodic_triangle Periodic_triangle;
typedef typename Base::Periodic_point Periodic_point;
typedef typename Base::Periodic_segment Periodic_segment;
typedef typename Base::Periodic_triangle Periodic_triangle;
typedef typename Base::Vertex Vertex;
typedef typename Base::Edge Edge;
typedef typename Base::Face Face;
typedef typename Base::Vertex Vertex;
typedef typename Base::Edge Edge;
typedef typename Base::Face Face;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Face_handle Face_handle;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Face_handle Face_handle;
typedef typename Base::size_type size_type;
typedef typename Base::difference_type difference_type;
typedef typename Base::size_type size_type;
typedef typename Base::difference_type difference_type;
typedef typename Base::Face_iterator Face_iterator;
typedef typename Base::Edge_iterator Edge_iterator;
typedef typename Base::Vertex_iterator Vertex_iterator;
typedef typename Base::Face_circulator Face_circulator;
typedef typename Base::Edge_circulator Edge_circulator;
typedef typename Base::Vertex_circulator Vertex_circulator;
typedef typename Base::Line_face_circulator Line_face_circulator;
typedef typename Base::Face_iterator Face_iterator;
typedef typename Base::Edge_iterator Edge_iterator;
typedef typename Base::Vertex_iterator Vertex_iterator;
typedef typename Base::Face_circulator Face_circulator;
typedef typename Base::Edge_circulator Edge_circulator;
typedef typename Base::Vertex_circulator Vertex_circulator;
typedef typename Base::Line_face_circulator Line_face_circulator;
private:
typedef typename GT::FT FT;
typedef typename GT::FT FT;
class Dummy_point {
private:
Point _pt;
bool _is_inserted;
Vertex_handle _vh;
class Dummy_point {
private:
Point _pt;
bool _is_inserted;
Vertex_handle _vh;
public:
public:
Dummy_point(FT x, FT y): _pt(x,y), _is_inserted(true) {}
Dummy_point(Point p): _pt(p), _is_inserted(true) {}
Dummy_point(FT x, FT y): _pt(x,y), _is_inserted(true) {}
Dummy_point(Point p): _pt(p), _is_inserted(true) {}
Point operator()() const { return _pt; }
bool is_inserted() const { return _is_inserted; }
Vertex_handle vertex() const { return _vh; }
void set_inserted(bool val) { _is_inserted = val; }
void set(Point val) { _pt = val; }
void set_vertex(Vertex_handle v) { _vh = v; }
};
Point operator()() const { return _pt; }
bool is_inserted() const { return _is_inserted; }
Vertex_handle vertex() const { return _vh; }
void set_inserted(bool val) { _is_inserted = val; }
void set(Point val) { _pt = val; }
void set_vertex(Vertex_handle v) { _vh = v; }
};
std::vector<Dummy_point> dummy_points;
std::vector<Dummy_point> dummy_points;
public:
typedef Point value_type;
typedef const value_type& const_reference;
typedef Tag_false Weighted_tag;
typedef Point value_type;
typedef const value_type& const_reference;
typedef Tag_false Weighted_tag;
protected:
int f_cnt, v_cnt, n_dpt;
int f_cnt, v_cnt, n_dpt;
private:
bool is_removable(Vertex_handle v, Delaunay_triangulation_2<GT, TDS>& dt, std::map<Vertex_handle, Vertex_handle>& vmap);
bool is_removable(Vertex_handle v, Delaunay_triangulation_2<GT, TDS>& dt, std::map<Vertex_handle, Vertex_handle>& vmap);
public:
Periodic_4_hyperbolic_Delaunay_triangulation_2(Geometric_traits gt) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(gt) { n_dpt = 14; }
Periodic_4_hyperbolic_Delaunay_triangulation_2(Geometric_traits gt) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(gt) { n_dpt = 14; }
Periodic_4_hyperbolic_Delaunay_triangulation_2(
const Circle_2 domain = Circle_2(Point_2(FT(0),FT(0)), FT(1*1)),
const Geometric_traits &gt = Geometric_traits() ) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(domain, gt) { n_dpt = 14; }
Periodic_4_hyperbolic_Delaunay_triangulation_2(
const Circle_2 domain = Circle_2(Point_2(FT(0),FT(0)), FT(1*1)),
const Geometric_traits &gt = Geometric_traits() ) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(domain, gt) { n_dpt = 14; }
Periodic_4_hyperbolic_Delaunay_triangulation_2(const Periodic_4_hyperbolic_Delaunay_triangulation_2& tr) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(tr) { }
Periodic_4_hyperbolic_Delaunay_triangulation_2(const Periodic_4_hyperbolic_Delaunay_triangulation_2& tr) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(tr) { }
std::vector<Vertex_handle> insert_dummy_points(bool rational = true);
std::vector<Vertex_handle> insert_dummy_points(bool rational = true);
Vertex_handle insert(const Point &p,
Face_handle start = Face_handle() );
Vertex_handle insert(const Point &p, Face_handle start = Face_handle() );
template < class InputIterator >
std::ptrdiff_t insert(InputIterator first,
InputIterator last);
template < class InputIterator >
std::ptrdiff_t
insert(InputIterator first, InputIterator last) {
size_type n = this->number_of_vertices();
Face_handle locate(const Point& p, Locate_type& lt, int& li, const Face_handle fh = Face_handle()) const {
Offset lo;
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
std::vector<Point> points(first, last);
spatial_sort(points.begin(), points.end(), geom_traits());
Face_handle f;
for (typename std::vector<Point>::const_iterator p = points.begin(), end = points.end(); p != end; ++p)
f = insert(*p, f)->face();
Face_handle locate(const Point& p, const Face_handle fh = Face_handle()) const {
Offset lo;
Locate_type lt;
int li;
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
return this->number_of_vertices() - n;
}
Face_handle periodic_locate(const Point& p, Locate_type& lt, int& li, Offset& lo, const Face_handle fh = Face_handle()) const {
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
Face_handle locate(const Point& p, Locate_type& lt, int& li, const Face_handle fh = Face_handle()) const {
Offset lo;
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
Point_2 get_dummy_point(int i) const {
return dummy_points[i]();
}
Face_handle locate(const Point& p, const Face_handle fh = Face_handle()) const {
Offset lo;
Locate_type lt;
int li;
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
void remove(Vertex_handle v);
Face_handle periodic_locate(const Point& p, Locate_type& lt, int& li, Offset& lo, const Face_handle fh = Face_handle()) const {
return this->hyperbolic_locate(p, lt, li, lo, fh);
}
int number_of_dummy_points() { return n_dpt; }
Point_2 get_dummy_point(int i) const {
return dummy_points[i]();
}
bool _side_of_octagon( const Face_handle& fh, const Offset& offset) const {
int cnt = 0;
typename GT::Side_of_fundamental_octagon side;
for (int j = 0; j < 3; j++) {
Offset o = offset.inverse().append(fh->vertex(j)->get_offset());
Point p = o.apply( fh->vertex(j)->point() );
if ( side(p) == CGAL::ON_UNBOUNDED_SIDE ) {
if ( p.y() + tan(CGAL_PI / FT(8))*p.x() > 0 ) {
cnt++;
} else {
}
}
}
return (cnt == 0);
}
void remove(Vertex_handle v);
int number_of_dummy_points() { return n_dpt; }
bool _side_of_octagon( const Face_handle& fh, const Offset& offset) const {
int cnt = 0;
typename GT::Side_of_fundamental_octagon side;
for (int j = 0; j < 3; j++) {
Offset o = offset.inverse().append(fh->vertex(j)->get_offset());
Point p = o.apply( fh->vertex(j)->point() );
if ( side(p) == CGAL::ON_UNBOUNDED_SIDE ) {
if ( p.y() + tan(CGAL_PI / FT(8))*p.x() > 0 ) {
cnt++;
} else {
}
}
}
return (cnt == 0);
}
void clear() {
Base::clear();
insert_dummy_points(true);
}
void clear() {
Base::clear();
insert_dummy_points(true);
}
}; // class Periodic_4_hyperbolic_Delaunay_triangulation_2
@ -229,8 +238,8 @@ hyperbolic_diameter(typename Gt::Circle_2 c) {
typedef std::pair<CGAL::Circular_arc_point_2<CircK>, unsigned> IsectOutput;
typedef CGAL::Dispatch_output_iterator<
CGAL::cpp11::tuple<IsectOutput>,
CGAL::cpp0x::tuple< std::back_insert_iterator<std::vector<IsectOutput> > > > Dispatcher;
CGAL::cpp11::tuple<IsectOutput>,
CGAL::cpp0x::tuple< std::back_insert_iterator<std::vector<IsectOutput> > > > Dispatcher;
std::vector<IsectOutput> res0, res1;
Dispatcher disp1 = CGAL::dispatch_output<IsectOutput>( std::back_inserter(res1) );
@ -245,11 +254,11 @@ hyperbolic_diameter(typename Gt::Circle_2 c) {
Circ2 c0(p0, 1);
if (ell.is_degenerate()) {
//cout << "\tThis is degenerate case!" << endl;
return 5.;
//cout << "\tThis is degenerate case!" << endl;
return 5.;
} else {
CGAL::intersection(c0, ell, disp0);
CGAL::intersection(c2, ell, disp1);
CGAL::intersection(c0, ell, disp0);
CGAL::intersection(c2, ell, disp1);
}
Point a(to_double(res0[0].first.x()), to_double(res0[0].first.y()));
Point b(to_double(res0[1].first.x()), to_double(res0[1].first.y()));
@ -288,14 +297,14 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2<Gt,Tds>& dt, std::map<Ver
std::vector<Vertex_handle> bdry_verts;
Face_circulator nbf(tds().incident_faces(v)), done(nbf);
do {
int idx = nbf->index(v);
Offset off = nbf->offset(idx).inverse();
off = off*nbf->offset(ccw(idx));
Vertex_handle thisv = nbf->vertex(ccw(idx));
bdry_verts.push_back(thisv);
Point pt = off.apply(thisv->point());
Vertex_handle new_v = dt.insert(pt);
vmap.insert(std::pair<Vertex_handle, Vertex_handle>(new_v, thisv));
int idx = nbf->index(v);
Offset off = nbf->offset(idx).inverse();
off = off*nbf->offset(ccw(idx));
Vertex_handle thisv = nbf->vertex(ccw(idx));
bdry_verts.push_back(thisv);
Point pt = off.apply(thisv->point());
Vertex_handle new_v = dt.insert(pt);
vmap.insert(std::pair<Vertex_handle, Vertex_handle>(new_v, thisv));
} while (++nbf != done);
@ -303,35 +312,35 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2<Gt,Tds>& dt, std::map<Ver
FT max_diam = FT(0);
for (Finite_Delaunay_faces_iterator fit = dt.finite_faces_begin(); fit != dt.finite_faces_end(); fit++) {
bool is_good = true;
for (int i = 0; i < 3; i++) {
Vertex_handle this_v = vmap[fit->vertex(i)];
Vertex_handle prev_v = bdry_verts[n_verts - 1];
Vertex_handle curr_v = bdry_verts[0];
for (int j = 1; curr_v != this_v; j = (j+1)%n_verts) {
prev_v = curr_v;
curr_v = bdry_verts[j];
}
if (vmap[fit->vertex(ccw(i))] == prev_v) {
is_good = false;
break;
}
}
if (is_good) {
Circle c(fit->vertex(0)->point(),
fit->vertex(1)->point(),
fit->vertex(2)->point());
FT diam = hyperbolic_diameter<Gt>(c);
if (max_diam < diam) {
max_diam = diam;
}
}
bool is_good = true;
for (int i = 0; i < 3; i++) {
Vertex_handle this_v = vmap[fit->vertex(i)];
Vertex_handle prev_v = bdry_verts[n_verts - 1];
Vertex_handle curr_v = bdry_verts[0];
for (int j = 1; curr_v != this_v; j = (j+1)%n_verts) {
prev_v = curr_v;
curr_v = bdry_verts[j];
}
if (vmap[fit->vertex(ccw(i))] == prev_v) {
is_good = false;
break;
}
}
if (is_good) {
Circle c(fit->vertex(0)->point(),
fit->vertex(1)->point(),
fit->vertex(2)->point());
FT diam = hyperbolic_diameter<Gt>(c);
if (max_diam < diam) {
max_diam = diam;
}
}
}
if (max_diam < lim) {
return true;
return true;
} else {
return false;
return false;
}
}
@ -340,70 +349,67 @@ template < class Gt, class Tds >
inline
typename Periodic_4_hyperbolic_Delaunay_triangulation_2<Gt, Tds>::Vertex_handle
Periodic_4_hyperbolic_Delaunay_triangulation_2<Gt, Tds>::
insert(const Point &p, Face_handle start) {
insert(const Point &p, Face_handle hint) {
Vertex_handle v;
Vertex_handle v;
typedef typename Gt::Side_of_fundamental_octagon Side_of_fundamental_octagon;
typedef typename Gt::Side_of_fundamental_octagon Side_of_fundamental_octagon;
Side_of_fundamental_octagon check = Side_of_fundamental_octagon();
CGAL::Bounded_side side = check(p);
Side_of_fundamental_octagon check = Side_of_fundamental_octagon();
CGAL::Bounded_side side = check(p);
if (side != CGAL::ON_UNBOUNDED_SIDE) {
Offset loff;
if ( start == Face_handle() ) {
Locate_type lt;
int li;
start = this->euclidean_locate(p, lt, li, loff);
if (lt == Periodic_4_hyperbolic_Delaunay_triangulation_2<Gt, Tds>::VERTEX) {
return Vertex_handle();
}
}
if (side != CGAL::ON_UNBOUNDED_SIDE) {
Offset loff;
Locate_type lt;
int li;
Face_handle start = this->euclidean_locate(p, lt, li, loff, hint);
if (lt == Periodic_4_hyperbolic_Delaunay_triangulation_2<Gt, Tds>::VERTEX) {
return Vertex_handle();
}
std::vector<Face_handle> faces;
this->find_conflicts(start, p, loff, std::back_inserter(faces));
v = this->insert_in_hole(p, faces.begin(), faces.end());
v->set_idx(v_cnt++);
std::vector<Face_handle> faces;
this->find_conflicts(start, p, loff, std::back_inserter(faces));
v = this->insert_in_hole(p, faces.begin(), faces.end());
v->set_idx(v_cnt++);
Face_circulator ifc = tds().incident_faces(v), done(ifc);
do {
ifc->restore_offsets(loff);
ifc->tds_data().clear();
ifc->make_canonical();
ifc->set_number(f_cnt++);
} while (++ifc != done);
Face_circulator ifc = tds().incident_faces(v), done(ifc);
do {
ifc->restore_offsets(loff);
ifc->tds_data().clear();
ifc->make_canonical();
ifc->set_number(f_cnt++);
} while (++ifc != done);
Vertex_circulator ivc = tds().incident_vertices(v), done_v(ivc);
do {
ivc->remove_offset();
} while (++ivc != done_v);
Vertex_circulator ivc = tds().incident_vertices(v), done_v(ivc);
do {
ivc->remove_offset();
} while (++ivc != done_v);
CGAL_triangulation_assertion(this->is_valid());
CGAL_triangulation_assertion(this->is_valid());
for (int i = 0; i < dummy_points.size(); i++) {
if (dummy_points[i].is_inserted()) {
typedef Delaunay_triangulation_2<Gt, Tds> Delaunay;
Delaunay dt;
std::map<Vertex_handle, Vertex_handle> vmap;
for (int i = 0; i < dummy_points.size(); i++) {
if (dummy_points[i].is_inserted()) {
typedef Delaunay_triangulation_2<Gt, Tds> Delaunay;
Delaunay dt;
std::map<Vertex_handle, Vertex_handle> vmap;
if (is_removable(dummy_points[i].vertex(), dt, vmap)) {
//cout << "Removing dummy point " << i << endl;
remove(dummy_points[i].vertex());
dummy_points[i].set_inserted(false);
}
}
}
if (is_removable(dummy_points[i].vertex(), dt, vmap)) {
remove(dummy_points[i].vertex());
dummy_points[i].set_inserted(false);
}
}
}
n_dpt = 0;
for (int i = 0; i < dummy_points.size(); i++) {
if (dummy_points[i].is_inserted())
n_dpt++;
}
n_dpt = 0;
for (int i = 0; i < dummy_points.size(); i++) {
if (dummy_points[i].is_inserted())
n_dpt++;
}
return v;
}
return v;
}
return Vertex_handle();
return Vertex_handle();
}
@ -420,126 +426,126 @@ remove(Vertex_handle v) {
if (is_removable(v, dt, vmap)) {
typedef typename Delaunay::Finite_faces_iterator Finite_Delaunay_faces_iterator;
typedef std::pair<Face_handle, int> Neighbor_pair;
typedef std::pair<Edge, Neighbor_pair> Edge_neighbor;
typedef typename Delaunay::Finite_faces_iterator Finite_Delaunay_faces_iterator;
typedef std::pair<Face_handle, int> Neighbor_pair;
typedef std::pair<Edge, Neighbor_pair> Edge_neighbor;
std::vector<Edge> bdry_edges;
std::vector<Vertex_handle> bdry_verts;
std::map<Edge, Neighbor_pair> bdry_nbrs;
std::vector<Edge> bdry_edges;
std::vector<Vertex_handle> bdry_verts;
std::map<Edge, Neighbor_pair> bdry_nbrs;
Face_circulator nb = tds().incident_faces(v), done(nb);
std::vector<Face_handle> nbrs;
do {
int idx = nb->index(v);
Edge e = Edge(nb, idx);
bdry_edges.push_back(e);
Face_handle nbf = nb->neighbor(idx);
int nidx = 0;
if (nbf->neighbor(1) == nb) nidx = 1;
if (nbf->neighbor(2) == nb) nidx = 2;
CGAL_triangulation_assertion(nbf->neighbor(nidx) == nb);
bdry_nbrs.insert(Edge_neighbor(e, Neighbor_pair(nbf, nidx)));
bdry_verts.push_back(nb->vertex(ccw(idx)));
Face_circulator nb = tds().incident_faces(v), done(nb);
std::vector<Face_handle> nbrs;
do {
int idx = nb->index(v);
Edge e = Edge(nb, idx);
bdry_edges.push_back(e);
Face_handle nbf = nb->neighbor(idx);
int nidx = 0;
if (nbf->neighbor(1) == nb) nidx = 1;
if (nbf->neighbor(2) == nb) nidx = 2;
CGAL_triangulation_assertion(nbf->neighbor(nidx) == nb);
bdry_nbrs.insert(Edge_neighbor(e, Neighbor_pair(nbf, nidx)));
bdry_verts.push_back(nb->vertex(ccw(idx)));
nb->store_offsets(nb->offset(idx).inverse());
nbrs.push_back(nb);
nb++;
} while(nb != done);
nb->store_offsets(nb->offset(idx).inverse());
nbrs.push_back(nb);
nb++;
} while(nb != done);
for (int i = 0; i < bdry_edges.size(); i++) {
Edge e = bdry_edges[i];
Face_handle f = e.first;
int j = e.second;
}
for (int i = 0; i < bdry_edges.size(); i++) {
Edge e = bdry_edges[i];
Face_handle f = e.first;
int j = e.second;
}
int n_verts = bdry_verts.size();
std::vector<Face_handle> new_f;
for (Finite_Delaunay_faces_iterator fit = dt.finite_faces_begin(); fit != dt.finite_faces_end(); fit++) {
bool is_good = true;
for (int i = 0; i < 3; i++) {
Vertex_handle this_v = vmap[fit->vertex(i)];
Vertex_handle prev_v = bdry_verts[n_verts - 1];
Vertex_handle curr_v = bdry_verts[0];
for (int j = 1; curr_v != this_v; j = (j+1)%n_verts) {
prev_v = curr_v;
curr_v = bdry_verts[j];
}
if (vmap[fit->vertex(ccw(i))] == prev_v) {
is_good = false;
break;
}
}
int n_verts = bdry_verts.size();
std::vector<Face_handle> new_f;
for (Finite_Delaunay_faces_iterator fit = dt.finite_faces_begin(); fit != dt.finite_faces_end(); fit++) {
bool is_good = true;
for (int i = 0; i < 3; i++) {
Vertex_handle this_v = vmap[fit->vertex(i)];
Vertex_handle prev_v = bdry_verts[n_verts - 1];
Vertex_handle curr_v = bdry_verts[0];
for (int j = 1; curr_v != this_v; j = (j+1)%n_verts) {
prev_v = curr_v;
curr_v = bdry_verts[j];
}
if (vmap[fit->vertex(ccw(i))] == prev_v) {
is_good = false;
break;
}
}
if (is_good) {
Face_handle f = tds().create_face();
f->set_number(f_cnt++);
for (int j = 0; j < 3; j++) {
f->set_vertex(j, vmap[fit->vertex(j)]);
}
new_f.push_back(f);
}
}
if (is_good) {
Face_handle f = tds().create_face();
f->set_number(f_cnt++);
for (int j = 0; j < 3; j++) {
f->set_vertex(j, vmap[fit->vertex(j)]);
}
new_f.push_back(f);
}
}
int internb = 0;
int bdrynb = 0;
for (int i = 0; i < new_f.size(); i++) {
for (int k = 0; k < 3; k++) {
bool found_bdry = false;
for (int j = 0; j < bdry_verts.size(); j++) {
if (new_f[i]->vertex(ccw(k)) == bdry_verts[j] &&
new_f[i]->vertex(cw(k)) == bdry_verts[(j+1)%n_verts]) {
found_bdry = true;
Neighbor_pair nb = bdry_nbrs[bdry_edges[j]];
Face_handle nbf = nb.first;
int nbidx = nb.second;
tds().set_adjacency(nbf, nbidx, new_f[i], k);
bdrynb++;
break;
}
}
if (!found_bdry) {
for (int l = 0; l < new_f.size(); l++) {
if (l == i) continue;
for (int j = 0; j < 3; j++) {
if (new_f[i]->vertex(ccw(k)) == new_f[l]->vertex(cw(j)) &&
new_f[i]->vertex(cw(k)) == new_f[l]->vertex(ccw(j)) ) {
tds().set_adjacency(new_f[i], k, new_f[l], j);
internb++;
break;
}
}
}
}
}
}
int internb = 0;
int bdrynb = 0;
for (int i = 0; i < new_f.size(); i++) {
for (int k = 0; k < 3; k++) {
bool found_bdry = false;
for (int j = 0; j < bdry_verts.size(); j++) {
if (new_f[i]->vertex(ccw(k)) == bdry_verts[j] &&
new_f[i]->vertex(cw(k)) == bdry_verts[(j+1)%n_verts]) {
found_bdry = true;
Neighbor_pair nb = bdry_nbrs[bdry_edges[j]];
Face_handle nbf = nb.first;
int nbidx = nb.second;
tds().set_adjacency(nbf, nbidx, new_f[i], k);
bdrynb++;
break;
}
}
if (!found_bdry) {
for (int l = 0; l < new_f.size(); l++) {
if (l == i) continue;
for (int j = 0; j < 3; j++) {
if (new_f[i]->vertex(ccw(k)) == new_f[l]->vertex(cw(j)) &&
new_f[i]->vertex(cw(k)) == new_f[l]->vertex(ccw(j)) ) {
tds().set_adjacency(new_f[i], k, new_f[l], j);
internb++;
break;
}
}
}
}
}
}
for (int j = 0; j < new_f.size(); j++) {
for (int i = 0; i < 3; i++) {
new_f[j]->vertex(i)->set_face(new_f[j]);
}
new_f[j]->restore_offsets();
new_f[j]->make_canonical();
}
for (int j = 0; j < new_f.size(); j++) {
for (int i = 0; i < 3; i++) {
new_f[j]->vertex(i)->set_face(new_f[j]);
}
new_f[j]->restore_offsets();
new_f[j]->make_canonical();
}
for (int j = 0; j < bdry_edges.size(); j++) {
Face_handle f = bdry_edges[j].first;
int i = bdry_edges[j].second;
f->vertex(ccw(i))->remove_offset();
}
for (int j = 0; j < bdry_edges.size(); j++) {
Face_handle f = bdry_edges[j].first;
int i = bdry_edges[j].second;
f->vertex(ccw(i))->remove_offset();
}
for (int i = 0; i < nbrs.size(); i++) {
tds().delete_face(nbrs[i]);
}
tds().delete_vertex(v);
for (int i = 0; i < nbrs.size(); i++) {
tds().delete_face(nbrs[i]);
}
tds().delete_vertex(v);
CGAL_triangulation_assertion(this->is_valid());
CGAL_triangulation_assertion(this->is_valid());
} else { // is not removable
cout << " -> vertex cannot be removed!" << endl;
cout << " -> vertex cannot be removed!" << endl;
}
}

View File

@ -885,13 +885,13 @@ template <class GT, class TDS>
typename TDS::Face_handle Periodic_4_hyperbolic_triangulation_2<GT, TDS>::
euclidean_locate(const Point& p, Locate_type& lt, int& li, Offset& loff, Face_handle f) const {
typedef typename GT::Side_of_fundamental_octagon Side_of_fundamental_octagon;
// typedef typename GT::Side_of_fundamental_octagon Side_of_fundamental_octagon;
Side_of_fundamental_octagon check = Side_of_fundamental_octagon();
CGAL::Bounded_side side = check(p);
if (side != ON_BOUNDED_SIDE) {
return Face_handle();
}
// Side_of_fundamental_octagon check = Side_of_fundamental_octagon();
// CGAL::Bounded_side side = check(p);
// if (side != ON_BOUNDED_SIDE) {
// return Face_handle();
// }
// Handle the case where an initial Face_handle is not given
if (f == Face_handle()) {

View File

@ -19,6 +19,7 @@ if ( CGAL_FOUND )
create_single_source_cgal_program( "test_p4ht2_dummy_points.cpp" )
create_single_source_cgal_program( "test_p4ht2_locate.cpp" )
create_single_source_cgal_program( "test_p4ht2_insertion.cpp" )
create_single_source_cgal_program( "test_p4ht2_insertion_euclidean.cpp" )
create_single_source_cgal_program( "test_p4ht2_removal.cpp" )
create_single_source_cgal_program( "test_p4ht2_remove_dummy_points.cpp" )
create_single_source_cgal_program( "test_p4ht2_complex.cpp" )

View File

@ -88,9 +88,10 @@ int main(int argc, char** argv) {
tr.insert_dummy_points(true);
CGAL::Timer t1;
t1.start();
for (int j = 0; j < pts.size(); j++) {
tr.insert(pts[j]);
}
tr.insert(pts.begin(), pts.end());
//for (int j = 0; j < pts.size(); j++) {
// tr.insert(pts[j]);
//}
t1.stop();
extime1 += t1.time();
cout << "DONE! (# of vertices = " << tr.number_of_vertices() << ", time = " << t1.time() << " secs)" << endl;
@ -109,9 +110,10 @@ int main(int argc, char** argv) {
Euclidean_triangulation etr;
CGAL::Timer t3;
t3.start();
for (int j = 0; j < pts.size(); j++) {
etr.insert(pts[j]);
}
etr.insert(pts.begin(), pts.end());
//for (int j = 0; j < pts.size(); j++) {
// etr.insert(pts[j]);
//}
t3.stop();
extime3 += t3.time();
cout << "DONE! (# of vertices = " << etr.number_of_vertices() << ", time = " << t3.time() << " secs)" << endl;

View File

@ -78,9 +78,10 @@ int main(int argc, char** argv) {
CGAL::Timer tt;
tt.start();
for (int j = 0; j < pts.size(); j++) {
tr.insert(pts[j]);
}
tr.insert(pts.begin(), pts.end());
//for (int j = 0; j < pts.size(); j++) {
// tr.insert(pts[j]);
//}
tt.stop();
cout << "DONE! (# of vertices = " << tr.number_of_vertices() << ", time = " << tt.time() << " secs)" << endl;
extime += tt.time();

View File

@ -0,0 +1,79 @@
#include <CGAL/basic.h>
#include <boost/tuple/tuple.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_smallint.hpp>
#include <boost/random/variate_generator.hpp>
#include <CGAL/point_generators_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
//#include <CGAL/Hyperbolic_random_points_in_disc_2.h>
//#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h>
//#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h>
//#include <CGAL/CORE_Expr.h>
//#include <CGAL/leda_real.h>
#include <CGAL/Cartesian.h>
#include <CGAL/determinant.h>
#include <CGAL/Timer.h>
//typedef CORE::Expr NT;
//typedef leda_real NT;
typedef double NT;
typedef CGAL::Cartesian<NT> Kernel;
//typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_traits_2<Kernel> Traits;
typedef CGAL::Delaunay_triangulation_2<Kernel> Triangulation;
//typedef Hyperbolic_octagon_translation_matrix<NT> Octagon_matrix;
typedef Kernel::Point_2 Point;
typedef Triangulation::Vertex_handle Vertex_handle;
//typedef Traits::Side_of_fundamental_octagon Side_of_fundamental_octagon;
typedef CGAL::Creator_uniform_2<double, Point> Creator;
using namespace std;
int main(int argc, char** argv) {
if (argc < 2) {
cout << "usage: " << argv[0] << " [number_of_points_to_insert] [optional: number_of_iterations]" << endl;
return -1;
}
int N = atoi(argv[1]);
int iters = 1;
if (argc == 3) {
iters = atoi(argv[2]);
}
cout << "---- for best results, make sure that you have compiled me in Release mode ----" << endl;
double extime = 0.0;
for (int exec = 1; exec <= iters; exec++) {
std::vector<Point> pts;
pts.reserve(N);
CGAL::Random_points_in_disc_2<Point, Creator> g(1.0);
CGAL::cpp11::copy_n( g, N, std::back_inserter(pts));
cout << "iteration " << exec << ": inserting into triangulation (rational dummy points)... "; cout.flush();
Triangulation tr;
//tr.insert_dummy_points(true);
CGAL::Timer tt;
tt.start();
tr.insert(pts.begin(), pts.end());
//for (int j = 0; j < pts.size(); j++) {
// tr.insert(pts[j]);
//}
tt.stop();
cout << "DONE! (# of vertices = " << tr.number_of_vertices() << ", time = " << tt.time() << " secs)" << endl;
extime += tt.time();
}
double avgtime = extime / (double)iters;
cout << "---------------------------------------" << endl;
cout << "Average execution time over " << iters << " iterations: " << avgtime << " secs" << endl << endl;
return 0;
}