songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

removal of dummy points implemented + test

This commit is contained in:
Iordan Iordanov 2016-12-13 15:50:13 +01:00
parent ef50f56dd0
commit 5ee11130b5
4 changed files with 233 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

143
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h
View File

@ -31,6 +31,16 @@
#include <CGAL/Hyperbolic_octagon_word_4.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Circular_kernel_intersections.h>
#include <iterator>
#include <CGAL/intersections.h>
#include <CGAL/result_of.h>
#include <CGAL/iterator.h>
#include <boost/bind.hpp>
namespace CGAL {
@ -113,6 +123,28 @@ namespace CGAL {
typedef Triple< Vertex_handle, Vertex_handle, Vertex_handle >
Vertex_triple;
class Dummy_point {
private:
Point _pt;
bool _is_inserted;
Vertex_handle _vh;
public:
Dummy_point(FT x, FT y): _pt(x,y), _is_inserted(true) {}
Dummy_point(Point p): _pt(p), _is_inserted(true) {}
Point operator()() { return _pt; }
bool is_inserted() { return _is_inserted; }
Vertex_handle vertex() { 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;
public:
typedef Periodic_4_hyperbolic_triangulation_triangle_iterator_2<Self>
Periodic_triangle_iterator;
@ -127,7 +159,7 @@ namespace CGAL {
protected:
mutable std::vector<Vertex_handle> v_offsets;
int f_cnt, v_cnt;
int f_cnt, v_cnt, n_dpt;
private:
@ -135,12 +167,12 @@ namespace CGAL {
public:
Periodic_4_hyperbolic_Delaunay_triangulation_2(Geometric_traits gt) :
Periodic_4_hyperbolic_triangulation_2<GT, TDS>(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) { }
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) { }
@ -156,6 +188,8 @@ namespace CGAL {
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;
@ -163,7 +197,7 @@ namespace CGAL {
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() + (CGAL_PI / FT(8))*p.x() > 0 ) {
if ( p.y() + tan(CGAL_PI / FT(8))*p.x() > 0 ) {
cnt++;
} else {
}
@ -175,6 +209,70 @@ namespace CGAL {
}; // class Periodic_4_hyperbolic_Delaunay_triangulation_2
template <class Gt>
typename Gt::FT dist(typename Gt::Point_2 p, typename Gt::Point_2 q) {
typename Gt::FT r = (p.x() - q.x())*(p.x() - q.x()) + (p.y() - q.y())*(p.y() - q.y());
return sqrt(r);
}
template <class Gt>
typename Gt::FT
hyperbolic_diameter(typename Gt::Circle_2 c) {
typedef typename Gt::FT FT;
typedef typename Gt::Kernel K;
typedef typename Gt::Point_2 Point;
typedef typename Gt::Line_2 Line;
typedef typename Gt::Circle_2 Circle;
typedef CGAL::Exact_circular_kernel_2 CircK;
typedef CGAL::Point_2<CircK> Pt2;
typedef CGAL::Circle_2<CircK> Circ2;
typedef CGAL::Line_2<CircK> Line2;
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;
std::vector<IsectOutput> res0, res1;
Dispatcher disp1 = CGAL::dispatch_output<IsectOutput>( std::back_inserter(res1) );
Dispatcher disp0 = CGAL::dispatch_output<IsectOutput>( std::back_inserter(res0) );
Pt2 ct(to_double(c.center().x()), to_double(c.center().y()));
Pt2 p0(0, 0);
double r = to_double(c.squared_radius());
Line2 ell(p0, ct);
Circ2 c2(ct, r);
Circ2 c0(p0, 1);
if (ell.is_degenerate()) {
//cout << "\tThis is degenerate case!" << endl;
return 5.;
} else {
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()));
Point p(to_double(res1[0].first.x()), to_double(res1[0].first.y()));
Point q(to_double(res1[1].first.x()), to_double(res1[1].first.y()));
FT aq = dist<Gt>(a, q);
FT pb = dist<Gt>(p, b);
FT ap = dist<Gt>(a, p);
FT qb = dist<Gt>(q, b);
//cout << "aq = " << aq << ", pb = " << pb << " | ap = " << ap << ", qb = " << qb << endl;
double hyperdist = fabs(log(to_double((aq*pb)/(ap*qb))));
return hyperdist;
}
template <class Gt, class Tds>
bool
@ -187,7 +285,10 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2<Gt,Tds>& dt, std::map<Ver
typedef typename Delaunay::Finite_faces_iterator Finite_Delaunay_faces_iterator;
// This is a rational approximation of the limit (1/4 of the squared systole)
FT lim = FT(584)/FT(1000);
Circle lc(dummy_points[0](),
dummy_points[13](),
Offset(1,6,3).apply(dummy_points[13]()));
FT lim = hyperbolic_diameter<Gt>(lc)/FT(2);
// This is the exact value of the limit (1/4 of the squared systole)
//FT lim = FT( pow(acosh(1. + sqrt(2.)),2.)/4. );
@ -207,7 +308,7 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2<Gt,Tds>& dt, std::map<Ver
int n_verts = bdry_verts.size();
FT max_sq_radius = FT(0);
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;
@ -228,14 +329,14 @@ is_removable(Vertex_handle v, Delaunay_triangulation_2<Gt,Tds>& dt, std::map<Ver
Circle c(fit->vertex(0)->point(),
fit->vertex(1)->point(),
fit->vertex(2)->point());
if (max_sq_radius < c.squared_radius()) {
max_sq_radius = c.squared_radius();
FT diam = hyperbolic_diameter<Gt>(c);
if (max_diam < diam) {
max_diam = diam;
}
}
}
if (max_sq_radius < lim) {
if (max_diam < lim) {
return true;
} else {
return false;
@ -287,6 +388,26 @@ insert(const Point &p, Face_handle start) {
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;
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);
}
}
}
n_dpt = 0;
for (int i = 0; i < dummy_points.size(); i++) {
if (dummy_points[i].is_inserted())
n_dpt++;
}
return v;
}
@ -426,7 +547,7 @@ remove(Vertex_handle v) {
CGAL_triangulation_assertion(this->is_valid());
} else { // is not removable
cout << "Vertex " << v->idx() << " cannot be removed!" << endl;
//cout << "Vertex " << v->idx() << " cannot be removed!" << endl;
}
}

53
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h
View File

@ -65,32 +65,30 @@ namespace CGAL {
f_cnt = fcount;
v_cnt = vcount;
std::vector<typename GT::Point_2> pts;
if (rational) {
// Push back the origin
pts.push_back(Point( FT(0), FT(0) ));
dummy_points.push_back(Dummy_point( FT(0), FT(0) ));
// Compute rational approximation of internal points
pts.push_back( Point( FT( 1)/FT(2), FT(-4)/FT(19) ) ); // 0
pts.push_back( Point( FT( 1)/FT(2), FT( 4)/FT(19) ) ); // 1
pts.push_back( Point( FT( 4)/FT(19), FT( 1)/FT(2) ) ); // 2
pts.push_back( Point( FT(-4)/FT(19), FT( 1)/FT(2) ) ); // 3
pts.push_back( Point( FT(-1)/FT(2), FT( 4)/FT(19) ) ); // 4
pts.push_back( Point( FT(-1)/FT(2), FT(-4)/FT(19) ) ); // 5
pts.push_back( Point( FT(-4)/FT(19), FT(-1)/FT(2) ) ); // 6
pts.push_back( Point( FT( 4)/FT(19), FT(-1)/FT(2) ) ); // 7
dummy_points.push_back( Dummy_point( FT( 1)/FT(2), FT(-4)/FT(19) ) ); // 0
dummy_points.push_back( Dummy_point( FT( 1)/FT(2), FT( 4)/FT(19) ) ); // 1
dummy_points.push_back( Dummy_point( FT( 4)/FT(19), FT( 1)/FT(2) ) ); // 2
dummy_points.push_back( Dummy_point( FT(-4)/FT(19), FT( 1)/FT(2) ) ); // 3
dummy_points.push_back( Dummy_point( FT(-1)/FT(2), FT( 4)/FT(19) ) ); // 4
dummy_points.push_back( Dummy_point( FT(-1)/FT(2), FT(-4)/FT(19) ) ); // 5
dummy_points.push_back( Dummy_point( FT(-4)/FT(19), FT(-1)/FT(2) ) ); // 6
dummy_points.push_back( Dummy_point( FT( 4)/FT(19), FT(-1)/FT(2) ) ); // 7
// Compute rational approximations of the midpoints
pts.push_back( Point( FT(-9)/FT(14), FT(0) ) ); // 4
pts.push_back( Point( FT(-5)/FT(11), FT(-5)/FT(11) ) ); // 5
pts.push_back( Point( FT(0), FT(-9)/FT(14) ) ); // 6
pts.push_back( Point( FT(5)/FT(11), FT(-5)/FT(11) ) ); // 7
dummy_points.push_back( Dummy_point( FT(-9)/FT(14), FT(0) ) ); // 4
dummy_points.push_back( Dummy_point( FT(-5)/FT(11), FT(-5)/FT(11) ) ); // 5
dummy_points.push_back( Dummy_point( FT(0), FT(-9)/FT(14) ) ); // 6
dummy_points.push_back( Dummy_point( FT(5)/FT(11), FT(-5)/FT(11) ) ); // 7
// The vertex v_0
pts.push_back( Point( FT(97)/FT(125), FT(-26)/FT(81) ) ); // 0
dummy_points.push_back( Dummy_point( FT(97)/FT(125), FT(-26)/FT(81) ) ); // 0
} else { // Algebraic dummy points
@ -99,7 +97,7 @@ namespace CGAL {
FT F1 = FT(1);
FT F2 = FT(2);
pts.push_back(Point( FT(0), FT(0) )); // origin
dummy_points.push_back(Dummy_point( FT(0), FT(0) )); // origin
FT i1(CGAL::sqrt(FT(2) - sqrt(FT(2))));
FT i2(FT(2)*i1);
@ -107,17 +105,17 @@ namespace CGAL {
FT i4(-1 + i3);
Point a0(FT( CGAL::sqrt( i4 )), FT(0)); // internal point
a0 = rot8(a0);
pts.push_back(a0);
dummy_points.push_back(Dummy_point(a0));
for (int k = 1; k < 8; k++) {
a0 = rotate(a0);
pts.push_back(a0);
dummy_points.push_back(Dummy_point(a0));
}
pts.push_back(Point(FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)), F0)); // μ(s_0)
pts.push_back(Point(FT( -CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2) ), FT(-CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_1)
pts.push_back(Point(F0, FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)))); // μ(s_2)
pts.push_back(Point(FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)), -FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_3)
pts.push_back(Point(FT( CGAL::sqrt(CGAL::sqrt(F2) + F1)/ F2 ), - FT( CGAL::sqrt(CGAL::sqrt(F2) - F1)/ F2) )); // v_0
dummy_points.push_back(Dummy_point(FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)), F0)); // μ(s_0)
dummy_points.push_back(Dummy_point(FT( -CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2) ), FT(-CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_1)
dummy_points.push_back(Dummy_point(F0, FT( -CGAL::sqrt(CGAL::sqrt(F2) - F1)))); // μ(s_2)
dummy_points.push_back(Dummy_point(FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)), -FT(CGAL::sqrt( (CGAL::sqrt(F2) - F1) / F2)) )); // μ(s_3)
dummy_points.push_back(Dummy_point(FT( CGAL::sqrt(CGAL::sqrt(F2) + F1)/ F2 ), - FT( CGAL::sqrt(CGAL::sqrt(F2) - F1)/ F2) )); // v_0
}
@ -199,7 +197,8 @@ namespace CGAL {
Vertex_handle vertices[14];
for (int i = 0; i < vcount; i++) {
vertices[i] = tds().create_vertex();
vertices[i]->set_point(pts[i]);
vertices[i]->set_point(dummy_points[i]());
dummy_points[i].set_vertex(vertices[i]);
vertices[i]->set_idx(i);
}
@ -212,11 +211,7 @@ namespace CGAL {
faces[i]->set_number(i);
}
cout << "Setting dimension!" << endl;
tds().set_dimension(2);
cout << "Dimension is now: " << tds().dimension() << endl;
for (int i = 0; i < 4; i++) {
int fn = i;

7
Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt
View File

@ -16,10 +16,11 @@ if ( CGAL_FOUND )
include_directories (BEFORE "../../include")
#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_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_removal.cpp" )
create_single_source_cgal_program( "test_p4ht2_remove_dummy_points.cpp" )
else()

73
Periodic_4_hyperbolic_triangulation_2/test/Periodic_4_hyperbolic_triangulation_2/test_p4ht2_remove_dummy_points.cpp Normal file
View File

@ -0,0 +1,73 @@
#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/Periodic_4_hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h>
#include <CGAL/Periodic_4_hyperbolic_triangulation_dummy_14.h>
#include <CGAL/CORE_Expr.h>
#include <CGAL/Cartesian.h>
#include <CGAL/determinant.h>
typedef CORE::Expr NT;
typedef CGAL::Cartesian<NT> Kernel;
typedef Kernel::FT FT;
typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_traits_2<Kernel> Traits;
typedef CGAL::Periodic_4_hyperbolic_Delaunay_triangulation_2<Traits> Triangulation;
typedef Triangulation::Face_iterator Face_iterator;
typedef Triangulation::Vertex_handle Vertex_handle;
typedef Triangulation::Point Point;
typedef Traits::Side_of_fundamental_octagon Side_of_fundamental_octagon;
typedef CGAL::Creator_uniform_2<NT, Point > Creator;
int main(int argc, char** argv) {
if (argc < 2) {
cout << "usage: " << argv[0] << " [number_of_iterations]" << endl;
return -1;
}
int iter = atoi(argv[1]);
CGAL::Random_points_in_disc_2<Point, Creator> g( 1.0 );
Side_of_fundamental_octagon pred;
int min = 1000;
int max = -1;
double mean = 0.0;
for (int j = 0; j < iter; j++) {
cout << "Iteration " << (j+1) << "/" << iter << "..." << endl;
Triangulation tr;
tr.insert_dummy_points();
assert(tr.is_valid(true));
int cnt = 0;
do {
Point pt = *g;
++g;
if (pred(pt) != CGAL::ON_UNBOUNDED_SIDE) {
tr.insert(pt);
cnt++;
}
} while (tr.number_of_dummy_points() > 0);
assert(tr.is_valid());
if (cnt > max)
max = cnt;
if (cnt < min)
min = cnt;
mean += cnt;
}
mean /= (double)iter;
cout << "Finished " << iter << " iterations!" << endl;
cout << "Minimum number of points inserted: " << min << endl;
cout << "Maximum number of points inserted: " << max << endl;
cout << "Average number of points inserted: " << mean << endl;
return 0;
}