diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/CMakeLists.txt b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/CMakeLists.txt index 31816a3b677..cb1cb5bc460 100644 --- a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/CMakeLists.txt +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/CMakeLists.txt @@ -19,7 +19,7 @@ include(${CGAL_USE_FILE}) find_package(Qt5 QUIET COMPONENTS Widgets) -include_directories (BEFORE ../../include include) +include_directories (BEFORE ../../include include ../../../Hyperbolic_triangulation_2/include) # ui files, created with Qt Designer qt5_wrap_ui( uis Periodic_4_hyperbolic_triangulation_2.ui ) @@ -30,9 +30,14 @@ qt5_add_resources ( RESOURCE_FILES resources/Periodic_4_hyperbolic_triangulation # cpp files -add_executable ( Periodic_4_Dirichlet_region_demo Periodic_4_Dirichlet_region_demo.cpp ) -qt5_use_modules( Periodic_4_Dirichlet_region_demo Widgets) -add_to_cached_list( CGAL_EXECUTABLE_TARGETS Periodic_4_Dirichlet_region_demo ) -target_link_libraries( Periodic_4_Dirichlet_region_demo ${CGAL_LIBRARIES} ${QT_LIBRARIES} ${CGAL_QT_LIBRARIES} ) +add_executable ( Periodic_4_hyperbolic_triangulation_2_demo Periodic_4_hyperbolic_triangulation_2_demo.cpp ) +qt5_use_modules( Periodic_4_hyperbolic_triangulation_2_demo Widgets) +add_to_cached_list( CGAL_EXECUTABLE_TARGETS Periodic_4_hyperbolic_triangulation_2_demo ) +target_link_libraries( Periodic_4_hyperbolic_triangulation_2_demo ${CGAL_LIBRARIES} ${QT_LIBRARIES} ${CGAL_QT_LIBRARIES} ) +add_executable ( Periodic_4_hyperbolic_Delaunay_triangulation_2_demo Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.cpp ) +qt5_use_modules( Periodic_4_hyperbolic_Delaunay_triangulation_2_demo Widgets) +add_to_cached_list( CGAL_EXECUTABLE_TARGETS Periodic_4_hyperbolic_Delaunay_triangulation_2_demo ) +target_link_libraries( Periodic_4_hyperbolic_Delaunay_triangulation_2_demo ${CGAL_LIBRARIES} ${QT_LIBRARIES} ${CGAL_QT_LIBRARIES} ) + diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.cpp b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.cpp new file mode 100644 index 00000000000..7201514d068 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.cpp @@ -0,0 +1,779 @@ +#include + +// CGAL headers +#include +#include +#include +#include + +#include + +// to be deleted +#include +// + +#include + +// Qt headers +#include +#include +#include +#include +#include +#include + +#include + +// GraphicsView items and event filters (input classes) +#include "CGAL/Qt/TriangulationCircumcircle.h" + +#include "CGAL/Qt/TriangulationMovingPoint.h" +#include "CGAL/Qt/TriangulationConflictZone.h" +#include "CGAL/Qt/TriangulationRemoveVertex.h" +#include "CGAL/Qt/TriangulationPointInputAndConflictZone.h" +//#include +#include + +// store color +#include +// visualise color +#include + +// unique words +//#include +#include +#include + +#include + +// dummy points +//#include + +// for filtering +#include + +// for viewportsBbox +#include + +// the two base classes +#include "ui_Delaunay_triangulation_2.h" +#include + +#define OPTION_INSERT_DUMMY_POINTS 0 +#define OPTION_INSERT_FIELD_DUMMY_POINTS 0 + + +typedef CGAL::Exact_predicates_inexact_constructions_kernel InR; +typedef CGAL::Exact_predicates_exact_constructions_kernel R; +typedef CGAL::Triangulation_hyperbolic_traits_2 K; +typedef K::Point_2 Point; +// keep color +typedef TranslationInfo Vb_info; +typedef CGAL::Triangulation_vertex_base_with_info_2< Vb_info, K > + Vb; +typedef CGAL::Triangulation_face_base_with_info_2 + Fb; +typedef CGAL::Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2 > + Delaunay; +typedef Delaunay::Vertex_handle Vertex_handle; + + + + + +struct PointsComparator { + static double eps; + + bool operator() (const Point& l, const Point& r) const + { + if(l.x() < r.x() - eps) { + return true; + } + if(l.x() < r.x() + eps) { + if(l.y() < r.y() - eps) { + return true; + } + } + return false; + } +}; + +double PointsComparator::eps = 0.0001; + +string glabels[] = { "a", "\\bar{b}", "c", "\\bar{d}", "\\bar{a}", "b", "\\bar{c}", "d" }; + + +void recurr(vector& v, vector g, int depth = 1) { + if (depth > 1) { + + recurr(v, g, depth-1); + + vector tmp; + vector tmpw; + for (int i = 0; i < v.size(); i++) { + tmp.push_back(v[i]); + } + + for (int i = 0; i < tmp.size(); i++) { + for (int j = 0; j < g.size(); j++) { + v.push_back(tmp[i]*g[j]); + } + } + } else { + for (int i = 0; i < g.size(); i++) { + v.push_back(g[i]); + } + } +} + + +void my_unique_words(std::vector& p, Point input, int length) { + std::vector g; + get_generators(g); + std::vector v; + recurr(v, g, length); + std::set s; + + for (int i = 0; i < v.size(); i++) { + s.insert( v[i] ); + } + + cout << "Original point and images: " << endl; + cout << input.x() << ", " << input.y() << endl; + for (set::iterator it = s.begin(); it != s.end(); it++) { + HyperbolicOctagonGroup m = *it; + pair res; + res = m.apply(to_double(input.x()), to_double(input.y())); + cout << res.first << ", " << res.second << endl; + p.push_back( Point(res.first, res.second) ); + } + +} + + + +void apply_unique_words(std::vector& points, Point input = Point(0, 0), double threshold = 10, int word_length = 6, double d = .998) +{ + + cout << "apply_unique_words called with threshold = " << threshold << ", word_length = " << word_length << ", d = " << d << endl; + + static vector unique_words; + static bool generated = false; + if(generated == false) { + generate_unique_words(unique_words, threshold, word_length); + generated = true; + } + + //points.resize(unique_words.size()); + pair res; + for(size_t i = 0; i < unique_words.size(); i++) { + pair res; + res = unique_words[i].apply(to_double(input.x()), to_double(input.y())); + + double dist = res.first*res.first + res.second*res.second; + if(dist < d) { + points.push_back( Point(res.first, res.second) ); + } + } +} + +void apply_unique_words_G(std::vector& points, Point input = Point(0, 0), double threshold = 6/*13.5*/, int word_length = 6/*20*/) +{ + static vector unique_words; + static bool generated = false; + if(generated == false) { + generate_unique_words(unique_words, threshold, word_length); + generated = true; + + // to generate all words + /* + ofstream fwords("m_w"); + fwords << "local words;\nPrint(\"words!\\n\");\nwords := [\n"; + for(size_t i = 0; i < unique_words.size() - 1; i++) { + fwords << unique_words[i].label << "," < indices; + while(findices >> index) { + indices.push_back(index); + } + cout << "indices size " << indices.size() << endl; + + pair res; + for(size_t i = 0; i < indices.size(); i++) { + pair res; + if(indices[i] < unique_words.size() /*&& unique_words[indices[i]].length() < 13.*/) { + res = unique_words[indices[i]].apply(to_double(input.x()), to_double(input.y())); + points.push_back(Point(res.first, res.second)); + } + } + cout << "nb of points to insert " << points.size() << endl; +} + + +class MainWindow : + public CGAL::Qt::DemosMainWindow, + public Ui::Delaunay_triangulation_2 +{ + Q_OBJECT + +private: + + int cidx; + std::vector ccol; + + Delaunay dt; + QGraphicsEllipseItem * disk; + QGraphicsScene scene; + + CGAL::Qt::TriangulationGraphicsItem * dgi; + CGAL::Qt::VoronoiGraphicsItem * vgi; + + // for drawing Voronoi diagram of the orbit of the origin + CGAL::Qt::VoronoiGraphicsItem * origin_vgi; + + CGAL::Qt::TriangulationMovingPoint * mp; + CGAL::Qt::TriangulationConflictZone * cz; + CGAL::Qt::TriangulationRemoveVertex * trv; + CGAL::Qt::TriangulationPointInputAndConflictZone * pi; + CGAL::Qt::TriangulationCircumcircle * tcc; +public: + MainWindow(); + +public slots: + + void processInput(CGAL::Object o); + + void on_actionMovingPoint_toggled(bool checked); + + void on_actionShowConflictZone_toggled(bool checked); + + void on_actionCircumcenter_toggled(bool checked); + + void on_actionShowDelaunay_toggled(bool checked); + + void on_actionShowVoronoi_toggled(bool checked); + + void on_actionInsertPoint_toggled(bool checked); + + void on_actionInsertRandomPoints_triggered(); + + void on_actionLoadPoints_triggered(); + + void on_actionSavePoints_triggered(); + + void on_actionClear_triggered(); + + void on_actionRecenter_triggered(); + + virtual void open(QString fileName); + +signals: + void changed(); +}; + + +MainWindow::MainWindow() + : DemosMainWindow(), dt(K(1)) +{ + + cidx = 0; + for (int i = 0; i < 10; i++) + ccol.push_back(i); + + setupUi(this); + + this->graphicsView->setAcceptDrops(false); + + // Add Poincaré disk + qreal origin_x = 0, origin_y = 0, radius = 1, diameter = 2*radius; + qreal left_top_corner_x = origin_x - radius; + qreal left_top_corner_y = origin_y - radius; + qreal width = diameter, height = diameter; + + // set background + qreal eps = 0.01; + QGraphicsRectItem* rect = new QGraphicsRectItem(left_top_corner_x - eps, left_top_corner_y - eps, width + 2*eps, height + 2*eps); + rect->setPen(Qt::NoPen); + rect->setBrush(Qt::white); + scene.addItem(rect); + + // set disk + disk = new QGraphicsEllipseItem(left_top_corner_x, left_top_corner_y, width, height); + QPen pen; // creates a default pen + pen.setWidth(0); + //pen.setBrush(Qt::black); + pen.setBrush(QColor(200, 200, 0)); + disk->setPen(pen); + + scene.addItem(disk); + + // another input point, instead of the origin + + double phi = CGAL_PI / 8.; + double psi = CGAL_PI / 3.; + double rho = std::sqrt(cos(psi)*cos(psi) - sin(phi)*sin(phi)); + + Point origin = Point(0, 0); + const Point a(cos(phi)*cos(phi + psi)/rho, sin(phi)*cos(phi + psi)/rho); + + + // dt to form the octagon tessellation + vector origin_orbit; + apply_unique_words(origin_orbit, origin, 10, 6); + origin_orbit.push_back(Point(0, 0)); + // for(long i = 0; i < origin_orbit.size(); i++) { + // cout << origin_orbit[i] << endl; + // } + cout << "nb of points on the orbit of the origin: " << origin_orbit.size() << endl; + + Delaunay* dtO = new Delaunay(K(1)); + dtO->insert(origin_orbit.begin(), origin_orbit.end()); + + origin_vgi = new CGAL::Qt::VoronoiGraphicsItem(dtO); + origin_vgi->setVisible(true); + + QObject::connect(this, SIGNAL(changed()), + origin_vgi, SLOT(modelChanged())); + + QColor br(149, 179, 179); + origin_vgi->setEdgesPen(QPen(br, 0, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin)); + scene.addItem(origin_vgi); + + + // Add a GraphicItem for the Delaunay triangulation + dgi = new CGAL::Qt::TriangulationGraphicsItem(&dt); + + QObject::connect(this, SIGNAL(changed()), + dgi, SLOT(modelChanged())); + + dgi->setVerticesPen(QPen(Qt::red, 3, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin)); + dgi->setEdgesPen(QPen(QColor(200, 200, 0), 0, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin)); + scene.addItem(dgi); + + // Add a GraphicItem for the Voronoi diagram + vgi = new CGAL::Qt::VoronoiGraphicsItem(&dt); + + QObject::connect(this, SIGNAL(changed()), + vgi, SLOT(modelChanged())); + + QColor brown(139, 69, 19); + vgi->setEdgesPen(QPen(brown, 0, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin)); + scene.addItem(vgi); + vgi->hide(); + + // Setup input handlers. They get events before the scene gets them + // and the input they generate is passed to the triangulation with + // the signal/slot mechanism + pi = new CGAL::Qt::TriangulationPointInputAndConflictZone(&scene, &dt, this ); + + QObject::connect(pi, SIGNAL(generate(CGAL::Object)), + this, SLOT(processInput(CGAL::Object))); + + mp = new CGAL::Qt::TriangulationMovingPoint(&dt, this); + // TriangulationMovingPoint emits a modelChanged() signal each + // time the moving point moves. + // The following connection is for the purpose of emitting changed(). + QObject::connect(mp, SIGNAL(modelChanged()), + this, SIGNAL(changed())); + + trv = new CGAL::Qt::TriangulationRemoveVertex(&dt, this); + QObject::connect(trv, SIGNAL(modelChanged()), + this, SIGNAL(changed())); + + tcc = new CGAL::Qt::TriangulationCircumcircle(&scene, &dt, this); + tcc->setPen(QPen(Qt::red, 0, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin)); + + cz = new CGAL::Qt::TriangulationConflictZone(&scene, &dt, this); + + // + // Manual handling of actions + // + + QObject::connect(this->actionQuit, SIGNAL(triggered()), + this, SLOT(close())); + + // We put mutually exclusive actions in an QActionGroup + QActionGroup* ag = new QActionGroup(this); + ag->addAction(this->actionInsertPoint); + ag->addAction(this->actionMovingPoint); + ag->addAction(this->actionCircumcenter); + ag->addAction(this->actionShowConflictZone); + + // Check two actions + this->actionInsertPoint->setChecked(true); + this->actionShowDelaunay->setChecked(true); + + // // + // // Setup the scene and the view + // // + scene.setItemIndexMethod(QGraphicsScene::NoIndex); + scene.setSceneRect(left_top_corner_x, left_top_corner_y, width, height); + this->graphicsView->setScene(&scene); + this->graphicsView->setMouseTracking(true); + this->graphicsView->shear(230, 230); + + // // The navigation adds zooming and translation functionality to the + // // QGraphicsView + this->addNavigation(this->graphicsView); + + this->setupStatusBar(); + this->setupOptionsMenu(); + this->addAboutDemo(":/cgal/help/about_Delaunay_triangulation_2.html"); + this->addAboutCGAL(); + + this->addRecentFiles(this->menuFile, this->actionQuit); + connect(this, SIGNAL(openRecentFile(QString)), + this, SLOT(open(QString))); + + +#if OPTION_INSERT_DUMMY_POINTS == 1 + + std::vector pts; + cout << "Inserting dummy points now! " << endl; + dt.insert_dummy_points(pts); + for (int i = 0; i < pts.size(); i++) { + processInput(make_object(pts[i])); + } + cout << "Dummy points inserted! " << endl; + emit(changed()); + +#endif + + +#if OPTION_INSERT_FIELD_DUMMY_POINTS == 1 + +Point mypt(1. - sqrt(2.), 1. - sqrt(2.)); +processInput(make_object(mypt)); +emit(changed()); + +#endif + +} + + +void +MainWindow::processInput(CGAL::Object o) +{ + + typedef CGAL::Exact_predicates_inexact_constructions_kernel GT; + typedef GT::Point_2 Point_2; + Point_2 pp; + if (CGAL::assign(pp, o)) { + Point tmp(pp.x(), pp.y()); + o = make_object(tmp); + } + + Point p; + if(CGAL::assign(p, o)){ + QPointF qp(CGAL::to_double(p.x()), CGAL::to_double(p.y())); + + // note that if the point is on the boundary then the disk contains the point + if(disk->contains(qp)){ + //dt.insert(p); + + //delete + vector points; + //apply_unique_words(points, p, 4, 1, .998); + my_unique_words(points, p, 1); + + points.push_back(p); + Vertex_handle v; + for(size_t j = 0; j < points.size() ; j++) { + v = dt.insert(points[j]); + v->info().setColor(ccol[cidx]); + } + cidx = (cidx + 1) % ccol.size(); + // + } + // delete + else { + static double phi = CGAL_PI / 8.; + static double psi = CGAL_PI / 3.; + static double rho = std::sqrt(cos(psi)*cos(psi) - sin(phi)*sin(phi)); + + static Point origin = Point(0, 0); + static Point a(cos(phi)*cos(phi + psi)/rho, sin(phi)*cos(phi + psi)/rho); + static Point b((cos(psi)-sin(phi))/rho, 0.); + + static Point current_point_a = origin; + static Point current_point_b = origin; + static double dT = 0.05; + + static double dx_a = dT * CGAL::to_double(a.x()); + static double dy_a = dT * CGAL::to_double(a.y()); + current_point_a = Point(current_point_a.x() + dx_a, current_point_a.y() + dy_a); + if(current_point_a.x() > a.x()) { + current_point_a = origin; + } + + static double dx_b = dT * CGAL::to_double(b.x()); + static double dy_b = dT * CGAL::to_double(b.y()); + current_point_b = Point(current_point_b.x() + dx_b, current_point_b.y() + dy_b); + std::cout << current_point_b.x() << " : " << current_point_b.y() << std::endl; + if(current_point_b.x() > b.x()) { + current_point_b = origin; + } + + vector points; + + // current_point_b <-> current_point_b + Point current_point = current_point_a; + apply_unique_words(points, current_point); + points.push_back(current_point); + dt.clear(); + dt.insert(points.begin(), points.end()); + + } + + } + emit(changed()); + + cout << "v = " << dt.number_of_vertices() << ", f = " << dt.number_of_faces() << endl; + +} + + + + + +/* + * Qt Automatic Connections + * http://doc.trolltech.com/4.4/designer-using-a-component.html#automatic-connections + * + * setupUi(this) generates connections to the slots named + * "on__" + */ +void +MainWindow::on_actionInsertPoint_toggled(bool checked) +{ + if(checked){ + scene.installEventFilter(pi); + scene.installEventFilter(trv); + } else { + scene.removeEventFilter(pi); + scene.removeEventFilter(trv); + } +} + + +void +MainWindow::on_actionMovingPoint_toggled(bool checked) +{ + + if(checked){ + scene.installEventFilter(mp); + } else { + scene.removeEventFilter(mp); + } +} + + +void +MainWindow::on_actionShowConflictZone_toggled(bool checked) +{ + + if(checked){ + scene.installEventFilter(cz); + } else { + scene.removeEventFilter(cz); + } +} + +void +MainWindow::on_actionCircumcenter_toggled(bool checked) +{ + if(checked){ + scene.installEventFilter(tcc); + tcc->show(); + } else { + scene.removeEventFilter(tcc); + tcc->hide(); + } +} + + +void +MainWindow::on_actionShowDelaunay_toggled(bool checked) +{ + dgi->setVisibleEdges(checked); +} + + +void +MainWindow::on_actionShowVoronoi_toggled(bool checked) +{ + vgi->setVisible(checked); +} + + +void +MainWindow::on_actionClear_triggered() +{ + dt.clear(); + emit(changed()); +} + + +void +MainWindow::on_actionInsertRandomPoints_triggered() +{ + QRectF rect = CGAL::Qt::viewportsBbox(&scene); + CGAL::Qt::Converter convert; + //Circle_2 isor = convert(rect); + //CGAL::Random_points_in_disc_2 pg(1.0); + bool ok = false; + const int number_of_points = + QInputDialog::getInt(this, + tr("Number of random points"), + tr("Enter number of random points"), + 100, + 0, + std::numeric_limits::max(), + 1, + &ok); + + if(!ok) { + return; + } + + // wait cursor + QApplication::setOverrideCursor(Qt::WaitCursor); + //std::vector points; + //points.reserve(number_of_points); + + typedef CGAL::Exact_predicates_inexact_constructions_kernel GT; + typedef GT::Point_2 Point_2; + typedef GT::FT FT; + + vector pts; + Hyperbolic_random_points_in_disc_2(pts, number_of_points); + + for(int i = 0; i < number_of_points; ++i){ + processInput(make_object(pts[i])); + //points.push_back(*pg++); + } + //dt.insert(points.begin(), points.end()); + // default cursor + QApplication::restoreOverrideCursor(); + emit(changed()); +} + + +void +MainWindow::on_actionLoadPoints_triggered() +{ + QString fileName = QFileDialog::getOpenFileName(this, + tr("Open Points file"), + "."); + if(! fileName.isEmpty()){ + open(fileName); + } +} + + +void +MainWindow::open(QString fileName) +{ + // wait cursor + QApplication::setOverrideCursor(Qt::WaitCursor); + std::ifstream ifs(qPrintable(fileName)); + + K::Point p; + std::vector points; + while(ifs >> p) { + points.push_back(p); + } + dt.insert(points.begin(), points.end()); + + // default cursor + QApplication::restoreOverrideCursor(); + this->addToRecentFiles(fileName); + actionRecenter->trigger(); + emit(changed()); + +} + +void +MainWindow::on_actionSavePoints_triggered() +{ +/* // dump points + QString fileName = QFileDialog::getSaveFileName(this, + tr("Save points"), + "."); + if(! fileName.isEmpty()){ + std::ofstream ofs(qPrintable(fileName)); + for(Delaunay::Finite_vertices_iterator + vit = dt.finite_vertices_begin(), + end = dt.finite_vertices_end(); + vit!= end; ++vit) + { + ofs << vit->point() << std::endl; + } + } +*/ + + // take a snapshot + std::cout << "snapshot..."; + + const QRect viewerRect = graphicsView->mapFromScene(scene.sceneRect()).boundingRect(); + const QRect imageRect = QRect(QPoint(0, 0), viewerRect.size()); + + QImage snapshot(imageRect.size(), QImage::Format_ARGB32);//QImage::Format_ARGB32_Premultiplied + + QPainter painter(&snapshot); + painter.setRenderHint(QPainter::Antialiasing); + //painter.setRenderHint(QPainter::SmoothPixmapTransform); + + graphicsView->render(&painter, imageRect, viewerRect); + bool saved = snapshot.save("mysnap.png", "PNG", 100); + assert(saved == true); + + std::cout << "done" << std::endl; +} + + +void +MainWindow::on_actionRecenter_triggered() +{ + this->graphicsView->setSceneRect(dgi->boundingRect()); + this->graphicsView->fitInView(dgi->boundingRect(), Qt::KeepAspectRatio); +} + + +#include "Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.moc" + +int main(int argc, char **argv) +{ + QApplication app(argc, argv); + + app.setOrganizationDomain("geometryfactory.com"); + app.setOrganizationName("GeometryFactory"); + app.setApplicationName("Delaunay_triangulation_2 demo"); + + // Import resources from libCGALQt4. + // See http://doc.trolltech.com/4.4/qdir.html#Q_INIT_RESOURCE + Q_INIT_RESOURCE(File); + Q_INIT_RESOURCE(Triangulation_2); + Q_INIT_RESOURCE(Input); + Q_INIT_RESOURCE(CGAL); + + MainWindow mainWindow; + mainWindow.show(); + + QStringList args = app.arguments(); + args.removeAt(0); + Q_FOREACH(QString filename, args) { + mainWindow.open(filename); + } + + return app.exec(); +} diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_Dirichlet_region_demo.cpp b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_triangulation_2_demo.cpp similarity index 90% rename from Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_Dirichlet_region_demo.cpp rename to Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_triangulation_2_demo.cpp index d78edf7c4e7..5170b58c3a5 100644 --- a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_Dirichlet_region_demo.cpp +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_triangulation_2_demo.cpp @@ -6,9 +6,7 @@ #include #include -/* #include */ - -#include +#include // to be deleted #include @@ -43,8 +41,8 @@ // unique words //#include -#include -#include +#include +#include #include @@ -67,24 +65,19 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel InR; typedef CGAL::Exact_predicates_exact_constructions_kernel R; typedef CGAL::Triangulation_hyperbolic_traits_2 K; - typedef K::Point_2 Point; // keep color typedef TranslationInfo Vb_info; - typedef CGAL::Triangulation_vertex_base_with_info_2< Vb_info, K > Vb; - typedef CGAL::Triangulation_face_base_with_info_2 Fb; - typedef CGAL::Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2 > -/* typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2 > */ Delaunay; - typedef Delaunay::Vertex_handle Vertex_handle; -//typedef CGAL::Delaunay_hyperbolic_triangulation_2 Delaunay; + + struct PointsComparator { static double eps; @@ -105,12 +98,64 @@ struct PointsComparator { double PointsComparator::eps = 0.0001; -void apply_unique_words(std::vector& points, Point input = Point(0, 0), double threshold = 20, int word_length = 6, double d = .999) +string glabels[] = { "a", "\\bar{b}", "c", "\\bar{d}", "\\bar{a}", "b", "\\bar{c}", "d" }; + + +void recurr(vector& v, vector g, int depth = 1) { + if (depth > 1) { + + recurr(v, g, depth-1); + + vector tmp; + vector tmpw; + for (int i = 0; i < v.size(); i++) { + tmp.push_back(v[i]); + } + + for (int i = 0; i < tmp.size(); i++) { + for (int j = 0; j < g.size(); j++) { + v.push_back(tmp[i]*g[j]); + } + } + } else { + for (int i = 0; i < g.size(); i++) { + v.push_back(g[i]); + } + } +} + + +void my_unique_words(std::vector& p, Point input, int length) { + std::vector g; + get_generators(g); + std::vector v; + recurr(v, g, length); + std::set s; + + for (int i = 0; i < v.size(); i++) { + s.insert( v[i] ); + } + + cout << "Original point and images: " << endl; + cout << input.x() << ", " << input.y() << endl; + for (set::iterator it = s.begin(); it != s.end(); it++) { + HyperbolicOctagonGroup m = *it; + pair res; + res = m.apply(to_double(input.x()), to_double(input.y())); + cout << res.first << ", " << res.second << endl; + p.push_back( Point(res.first, res.second) ); + } + +} + + + +void apply_unique_words(std::vector& points, Point input = Point(0, 0), double threshold = 10, int word_length = 6, double d = .998) { cout << "apply_unique_words called with threshold = " << threshold << ", word_length = " << word_length << ", d = " << d << endl; - static vector unique_words; + static vector unique_words; static bool generated = false; if(generated == false) { generate_unique_words(unique_words, threshold, word_length); @@ -132,7 +177,7 @@ void apply_unique_words(std::vector& points, Point input = Point(0, 0), d void apply_unique_words_G(std::vector& points, Point input = Point(0, 0), double threshold = 6/*13.5*/, int word_length = 6/*20*/) { - static vector unique_words; + static vector unique_words; static bool generated = false; if(generated == false) { generate_unique_words(unique_words, threshold, word_length); @@ -281,7 +326,7 @@ MainWindow::MainWindow() // dt to form the octagon tessellation vector origin_orbit; - apply_unique_words(origin_orbit, origin, 20, 6); + apply_unique_words(origin_orbit, origin, 10, 6); origin_orbit.push_back(Point(0, 0)); // for(long i = 0; i < origin_orbit.size(); i++) { // cout << origin_orbit[i] << endl; @@ -372,36 +417,11 @@ MainWindow::MainWindow() scene.setSceneRect(left_top_corner_x, left_top_corner_y, width, height); this->graphicsView->setScene(&scene); this->graphicsView->setMouseTracking(true); - - // // we want to adjust the coordinates of QGraphicsView to the coordinates of QGraphicsScene - // // the following line must do this: - //this->graphicsView->fitInView( scene.sceneRect(), Qt::KeepAspectRatio); - // // It does not do this sufficiently well. - // // Current solution: - - // QRect viewerRect = graphicsView->mapFromScene(scene.sceneRect()).boundingRect(); - // std::cout << "resolution before " << viewerRect.width() << " " << viewerRect.height() << std::endl; - - // // shear 230 230 - //this->graphicsView->shear(407, 407); - //this->graphicsView->shear(10, 10); - - // viewerRect = graphicsView->mapFromSthis->graphicsView->rotate(90);cene(scene.sceneRect()).boundingRect(); - // std::cout << "resolution after " << viewerRect.width() << " " << viewerRect.height() << std::endl; - - // // Turn the vertical axis upside down - //this->graphicsView->matrix().scale(1, -1); - - //this->graphicsView->scale(10, 10); - //this->graphicsView->centerOn(0.0, 0.0); - //this->graphicsView->translate(1.0, -20.0); this->graphicsView->shear(230, 230); - //this->graphicsView->rotate(90); // // The navigation adds zooming and translation functionality to the // // QGraphicsView this->addNavigation(this->graphicsView); - this->setupStatusBar(); this->setupOptionsMenu(); this->addAboutDemo(":/cgal/help/about_Delaunay_triangulation_2.html"); @@ -450,10 +470,12 @@ MainWindow::processInput(CGAL::Object o) //delete vector points; - apply_unique_words(points, p); + //apply_unique_words(points, p, 4, 1, .998); + my_unique_words(points, p, 1); + points.push_back(p); Vertex_handle v; - for(size_t j = 0; j < points.size(); j++) { + for(size_t j = 0; j < points.size() ; j++) { v = dt.insert(points[j]); v->info().setColor(ccol[cidx]); } @@ -715,7 +737,7 @@ MainWindow::on_actionRecenter_triggered() } -#include "Periodic_4_Dirichlet_region_demo.moc" +#include "Periodic_4_hyperbolic_triangulation_2_demo.moc" int main(int argc, char **argv) { diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt index 3ac6c2c4891..09983354192 100644 --- a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt @@ -20,7 +20,7 @@ include(${CGAL_USE_FILE}) find_package(Qt5 QUIET COMPONENTS Widgets) -include_directories (BEFORE ../../include include) +include_directories (BEFORE ../../include include ../../../Hyperbolic_triangulation_2/include) # ui files, created with Qt Designer qt5_wrap_ui( uis Periodic_4_hyperbolic_triangulation_2.ui ) diff --git a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/Periodic_4_hyperbolic_triangulation_2_demo.cpp b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/Periodic_4_hyperbolic_triangulation_2_demo.cpp index b4e2538fd0c..67783e71b04 100644 --- a/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/Periodic_4_hyperbolic_triangulation_2_demo.cpp +++ b/Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/Periodic_4_hyperbolic_triangulation_2_demo.cpp @@ -3,8 +3,8 @@ // CGAL headers #include #include -//#include -#include +#include +//#include #include // to be deleted @@ -43,7 +43,8 @@ typedef CGAL::Triangulation_hyperbolic_traits_2 K; typedef K::Point_2 Point_2; typedef K::Iso_rectangle_2 Iso_rectangle_2; -typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2 Delaunay; +//typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2 Delaunay; +typedef CGAL::Delaunay_hyperbolic_triangulation_2 Delaunay; class MainWindow : public CGAL::Qt::DemosMainWindow, @@ -332,7 +333,7 @@ MainWindow::on_actionModify_recursion_depth_triggered() return; } - dt.Set_recursion_depth(depth); + //dt.Set_recursion_depth(depth); } diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Base_hyperbolic_translations.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Base_hyperbolic_translations.h deleted file mode 100644 index d31a8468a30..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Base_hyperbolic_translations.h +++ /dev/null @@ -1,58 +0,0 @@ -#ifndef CGAL_HYPERBOLIC_BASE_TRANSLATIONS_2_H -#define CGAL_HYPERBOLIC_BASE_TRANSLATIONS_2_H - -#include - -template -struct Element -{ - typedef typename std::list List; - typedef CGAL::Circulator_from_container Circulator; - - Translation g; - - // circulator iterator to an inverse translation in the list - Circulator inverse; -}; - -template -class Base_hyperbolic_translations -{ -public: - typedef typename Gt::FT FT; - typedef typename std::complex complex; - typedef CGAL::Hyperbolic_isometry_2 Hyperbolic_isometry; - - //typedef Element Element_t; - //typedef std::list List; - //typedef typename List::iterator List_iterator; - //typedef CGAL::Circulator_from_container Circulator; - - typedef std::pair Node; - typedef std::vector Vector; - typedef typename Vector::iterator Vector_iterator; - - - virtual Hyperbolic_isometry& a() = 0; - virtual Hyperbolic_isometry& b() = 0; - virtual Hyperbolic_isometry& c() = 0; - virtual Hyperbolic_isometry& d() = 0; - - const Vector& get_vector_of_translations(); - Vector_iterator vector_begin(); - Vector_iterator vector_end(); - //List_iterator list_begin(); - //List_iterator list_end(); - //List& list(); - -private: - - void compute_g(); - //void compute_l(); - void compute(); - - Vector g; - //List l; -}; - -#endif // CGAL_HYPERBOLIC_BASE_TRANSLATIONS_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Cross_translations.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Cross_translations.h deleted file mode 100644 index af6c5a3550b..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Cross_translations.h +++ /dev/null @@ -1,125 +0,0 @@ -#ifndef CGAL_HYPERBOLIC_CROSS_TRANSLATIONS_2_H -#define CGAL_HYPERBOLIC_CROSS_TRANSLATIONS_2_H - -#include - - -template -class Cross_translations -{ -public: - - enum Direction { - A = 0, // 0 - B, // 1 - InvA, // 2 - InvB, // 3 - C, // 4 - D, // 5 - InvC, // 6 - InvD // 7 - }; - - typedef typename Gt::FT FT; - typedef typename std::complex complex; - typedef CGAL::Hyperbolic_isometry_2 Hyperbolic_isometry; - typedef std::vector Vector; - typedef typename Vector::iterator Vector_iterator; - - - Cross_translations() { - computed = false; - compute(); - } - - Hyperbolic_isometry& a() - { - compute(); - return g[A]; - } - - Hyperbolic_isometry& b() - { - compute(); - return g[B]; - } - - Hyperbolic_isometry& c() - { - compute(); - return g[C]; - } - - Hyperbolic_isometry& d() - { - compute(); - return g[D]; - } - - const Vector& get_vector_of_translations() - { - compute(); - return g; - } - - Vector_iterator vector_begin() - { - compute(); - return g.begin(); - } - - Vector_iterator vector_end() - { - compute(); - return g.end(); - } - - -private: - - void compute_g() - { - - const FT k1 = (FT(2) + CGAL::sqrt(2.))/FT(2); - const FT k2 = CGAL::sqrt(CGAL::sqrt(2.)); - const FT k3 = (CGAL::sqrt(2.)*k2)/FT(2); - - std::complex m(k1, k1); - std::complex n(k2*k1, -k3); - - g.resize(8); - - // a - g[A] = Hyperbolic_isometry(conj(m), conj(n)); - g[InvA] = g[A].inverse(); - - // b - g[B] = Hyperbolic_isometry(m, -n); - g[InvB] = g[B].inverse(); - - // c - g[C] = Hyperbolic_isometry(conj(m), -conj(n)); - g[InvC] = g[C].inverse(); - - // d - g[D] = Hyperbolic_isometry(m, n); - g[InvD] = g[D].inverse(); - - } - - - void compute() - { - if(!computed) { - compute_g(); - computed = true; - } - } - - bool computed; - Vector g; - -}; - - -#endif // CGAL_HYPERBOLIC_CROSS_TRANSLATIONS_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h index 902625a1349..1e6b53060e2 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h @@ -134,7 +134,6 @@ public: : Delaunay_triangulation_2(tr) { CGAL_triangulation_postcondition( this->is_valid() );} - void mark_star(Vertex_handle v) const { if(!is_star_bounded(v)) { diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Diametric_translations.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Diametric_translations.h deleted file mode 100644 index 4dc5065759d..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Diametric_translations.h +++ /dev/null @@ -1,132 +0,0 @@ -#ifndef CGAL_HYPERBOLIC_DIAMETRIC_TRANSLATIONS_2_H -#define CGAL_HYPERBOLIC_DIAMETRIC_TRANSLATIONS_2_H - -#include - -template -class Diametric_translations -{ -public: - - enum Direction { - A = 0, // 0 - InvB, // 1 - C, // 2 - InvD, // 3 - InvA, // 4 - B, // 5 - InvC, // 6 - D // 7 - }; - - typedef typename Gt::FT FT; - typedef typename std::complex complex; - typedef CGAL::Hyperbolic_isometry_2 Hyperbolic_isometry; - typedef std::vector Vector; - typedef typename Vector::iterator Vector_iterator; - - - Diametric_translations() { - computed = false; - compute(); - } - - Hyperbolic_isometry& a() - { - compute(); - return g[A]; - } - - Hyperbolic_isometry& b() - { - compute(); - return g[B]; - } - - Hyperbolic_isometry& c() - { - compute(); - return g[C]; - } - - Hyperbolic_isometry& d() - { - compute(); - return g[D]; - } - - const Vector& get_vector_of_translations() - { - compute(); - return g; - } - - Vector_iterator vector_begin() - { - compute(); - return g.begin(); - } - - Vector_iterator vector_end() - { - compute(); - return g.end(); - } - - -private: - - void compute_g() - { - - const FT k1 = FT(1) + CGAL::sqrt(2.); - const FT k2 = FT(0); - - std::complex m(k1, k2); - - // This is the multiplicative factor for all b's - const FT k3 = CGAL::sqrt(2.) * CGAL::sqrt( k1 ); - - std::vector< std::complex > n; - - // Euler's identity: exp(i k \theta) = cos(k \theta) + i sin(k \theta) - for (int kk = 0; kk < 8; kk++) { - n.push_back( std::complex( k3 * cos( FT(kk) * CGAL_PI / FT(4)), k3 * sin( FT(kk) * CGAL_PI / FT(4) ) ) ); - } - - g.resize(8); - - // a - g[A] = Hyperbolic_isometry(m, n[A]); - g[InvA] = g[A].inverse(); - - // b - g[B] = Hyperbolic_isometry(m, n[B]); - g[InvB] = g[B].inverse(); - - // c - g[C] = Hyperbolic_isometry(m, n[C]); - g[InvC] = g[C].inverse(); - - // d - g[D] = Hyperbolic_isometry(m, n[D]); - g[InvD] = g[D].inverse(); - - } - - - void compute() - { - if(!computed) { - compute_g(); - computed = true; - } - } - - bool computed; - Vector g; - -}; - - -#endif // CGAL_HYPERBOLIC_DIAMETRIC_TRANSLATIONS_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h deleted file mode 100644 index e69de29bb2d..00000000000 diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_isometry_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_isometry_2.h deleted file mode 100644 index 1d420f26a79..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_isometry_2.h +++ /dev/null @@ -1,105 +0,0 @@ -// Copyright (c) 2011 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org); you may redistribute it under -// the terms of the Q Public License version 1.0. -// See the file LICENSE.QPL distributed with CGAL. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL: svn+ssh://mbogdanov@scm.gforge.inria.fr/svn/cgal/trunk/Triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h $ -// $Id: Hyperboloic_isometry_2.h 57323 2010-07-05 10:07:39Z sloriot $ -// -// -// Author(s) : Mikhail Bogdanov - -#ifndef CGAL_HYPERBOLIC_ISOMETRY_2_H -#define CGAL_HYPERBOLIC_ISOMETRY_2_H - -#include - -#include - -namespace CGAL { - -/* - z -> (m*z + n)/(conj(n)*z + conj(m)) - */ -template -class Hyperbolic_isometry_2 -{ -public: - typedef Gt Geom_traits; - typedef typename Gt::FT FT; - typedef typename Gt::Point Point; - typedef std::complex complex; - - Hyperbolic_isometry_2(const complex& m = complex(1, 0), const complex& n = complex(0, 0)) - : m_(m), n_(n) - { - //std::cout << "Creating hyperbolic isometry with m = " << m << " and n = " << n << std::endl; - } - - ~Hyperbolic_isometry_2() - { - } - - complex DoAction(const complex& z) const - { - //std::cout << "Hyperbolic isometry applied!" << std::endl; - return (m_*z + n_)/(std::conj(n_)*z + std::conj(m_)); - } - - Point DoAction(const Point& p) const - { - complex z = complex(p.x(), p.y()); - complex result = DoAction(z); - - return Point(result.real(), result.imag()); - } - - Hyperbolic_isometry_2 inverse() const - { - return Hyperbolic_isometry_2(-std::conj(m_), n_); - } - - Hyperbolic_isometry_2 operator * (const Hyperbolic_isometry_2& other) const - { - return Hyperbolic_isometry_2(m_*other.m_ + n_*std::conj(other.n_), m_*other.n_ + n_*std::conj(other.m_)); - } - - void file_output(std::ostream& os) const - { - os << m_ << " " << n_ << std::endl; - } - - complex m() const - { - return m_; - } - - complex n() const - { - return n_; - } - -private: - complex m_; - complex n_; -}; - -template -std::ostream& -operator<<(std::ostream& os, const Hyperbolic_isometry_2 &isometry) -{ - isometry.file_output(os); - return os; -} - -} // namespace CGAL - -#endif // CGAL_HYPERBOLIC_ISOMETRY_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_octagon_group.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_octagon_group.h new file mode 100644 index 00000000000..238fc726390 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_octagon_group.h @@ -0,0 +1,575 @@ +// Copyright (c) 2010 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// +// Author(s) : Mikhail Bogdanov + +#ifndef CGAL_HYPERBOLIC_OCTAGON_GROUP_H +#define CGAL_HYPERBOLIC_OCTAGON_GROUP_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class Hyperbolic_octagon_group +{ + +public: + static Square_root_2_field factor; // The multiplicative factor present in the general formula: sqrt(2) - 1 + + typedef Hyperbolic_octagon_group Self; + typedef complex Matrix_element; + + Matrix_element A; + Matrix_element B; + string label; + + Hyperbolic_octagon_group(const Matrix_element& A_, const Matrix_element& B_, + const string& label_ = string("") ) : + A(A_), B(B_), label(label_) {} + + string merge_labels(const Self& rh) const + { + return label + "*" + rh.label; + } + + Self operator*(const Self& rh) const + { + return Self( A*rh.A + factor*B*conj(rh.B), + A*rh.B + B*conj(rh.A), merge_labels(rh) ); + } + + Self inverse() const + { + Self inv = Self(conj(A), -B); + string inv_label; + for(long i = label.size() - 1; i >= 0; i--) { + if(label[i] >= 'a' && label[i] <= 'd') { + inv_label.push_back(label[i]); + inv_label.push_back('^'); + inv_label.push_back('-'); + inv_label.push_back('1'); + } + if(label[i] == '*') { + inv_label.push_back('*'); + } + if(label[i] == '1') { + assert(i - 3 >= 0); + assert(label[i - 3] >= 'a' && label[i - 3] <= 'd'); + inv_label.push_back(label[i - 3]); + i = i - 3; + } + } + inv.label = inv_label; + return inv; + } + + // rotation \pi/4 + Hyperbolic_octagon_group rotate() const + { + Square_root_2_field B1 = real(B); + Square_root_2_field B2 = imag(B); + + // sqrt(2) + Square_root_2_field k = Square_root_2_field(0, 1); + + Square_root_2_field BB1 = (B1 - B2)*k; + Square_root_2_field BB2 = (B1 + B2)*k; + + assert(BB2.l % 2 == 0 && BB2.r % 2 == 0); + + BB1.l = BB1.l/2; + BB1.r = BB1.r/2; + BB2.l = BB2.l/2; + BB2.r = BB2.r/2; + + return Hyperbolic_octagon_group(A, Matrix_element(BB1, BB2)); + } + + Square_root_2_field trace() const + { + return Square_root_2_field(2, 0)*real(A); + } + + double length() const + { + typedef long double ld; + + ld l = real(A).l; + ld r = real(A).r; + ld tr = l + sqrt(2.)*r; + if (tr < 0) { + tr = -tr; + } + + return 2.*acosh(tr); + } + + // determinant == 1 + Matrix_element det() const + { + return norm(A) - factor * norm(B); + } + + static complex toComplexDouble(Matrix_element M) //const + { + Square_root_2_field rl = real(M); + Square_root_2_field img = imag(M); + + return complex(rl.l + sqrt(2.)*rl.r, img.l + sqrt(2.)*img.r); + } + + pair apply(double x, double y) + { + typedef complex Cmpl; + Cmpl Aa = toComplexDouble(A); + Cmpl Bb = toComplexDouble(B); + + double ax = sqrt(factor.l + sqrt(2.)*factor.r); + + Cmpl z(x, y); + Cmpl res = (Aa*z + ax*Bb)/(ax*(conj(Bb)*z) + conj(Aa)); + return pair(real(res), imag(res)); + } + +}; + + +template +Square_root_2_field Hyperbolic_octagon_group::factor = Square_root_2_field(-1, 1); + +// just to give an order(ing) +template +bool operator < (const complex >& lh, + const complex >& rh) +{ + if (real(lh) < real(rh)) { + return true; + } + + if (real(lh) == real(rh)) { + if (imag(lh) < imag(rh)) { + return true; + } + } + + return false; +} + +// just to order octagon_matrices +template +bool operator < (const Hyperbolic_octagon_group& lh, + const Hyperbolic_octagon_group& rh) +{ + if (lh.A < rh.A) { + return true; + } + + if (lh.A == rh.A ) { + if (lh.B < rh.B) { + return true; + } + } + + return false; +} + +template +bool operator == (const Hyperbolic_octagon_group& lh, + const Hyperbolic_octagon_group& rh) +{ + return (lh.A == rh.A && lh.B == rh.B); +} + +template +ostream& operator<<(ostream& os, const Hyperbolic_octagon_group& m) +{ + os << m.A << " " << m.B; + return os; +} + +typedef long long ll; +typedef Square_root_2_field SqrtField; +typedef Hyperbolic_octagon_group HyperbolicOctagonGroup; +typedef HyperbolicOctagonGroup::Matrix_element Element; + + +enum Direction { + G_A = 0, // 0 + G_InvB, // 1 + G_C, // 2 + G_InvD, // 3 + G_InvA, // 4 + G_B, // 5 + G_InvC, // 6 + G_D // 7 +}; + + +void get_generators(vector& gens) +{ + // This is a in the matrix, equal to sqrt(2) + 1 + Element A = Element(SqrtField(1, 1), SqrtField(0, 0)); + + // This vector holds all other elements, results of the exponentials for various k + vector B(8, Element(SqrtField(0, 0), SqrtField(0, 0))); + + // Corrected ordering (initial ordering is present in backup file) + B[G_A] = A * Element(SqrtField( 0, 0), SqrtField( 0, -1)); + B[G_B] = A * Element(SqrtField( 1, 0), SqrtField(-1, 0)); + B[G_C] = A * Element(SqrtField( 0, 1), SqrtField( 0, 0)); + B[G_D] = A * Element(SqrtField( 1, 0), SqrtField( 1, 0)); + B[G_InvA] = A * Element(SqrtField( 0, 0), SqrtField( 0, 1)); + B[G_InvB] = A * Element(SqrtField(-1, 0), SqrtField( 1, 0)); + B[G_InvC] = A * Element(SqrtField( 0, -1), SqrtField( 0, 0)); + B[G_InvD] = A * Element(SqrtField(-1, 0), SqrtField(-1, 0)); + + + + string labels[8] = {string("a"), string("\\bar{b}"), string("c"), string("\\bar{d}"), + string("\\bar{a}"), string("b"), string("\\bar{c}"), string("d")}; + for(int i = 0; i < 8; i++) { + gens.push_back(HyperbolicOctagonGroup(A, B[i], labels[i])); + } +} + +// a, \bar{b}, c, \bar{d}, \bar{a}, b, \bar{c}, d +vector gens; + +bool IsCanonical(const HyperbolicOctagonGroup& m); + +void generate_words( set& words, vector& prev, int depth, double threshold ) +{ + if (depth == 1) { + for(int i = 0; i < 8; i++) { + words.insert( gens[i] ); + prev.push_back( gens[i] ); + } + return; + } + + vector els; + generate_words( words, els, depth - 1, threshold); + + HyperbolicOctagonGroup temp = HyperbolicOctagonGroup(Element(), Element()); + ll size = els.size(); + bool is_new = false; + for(ll k = 0; k < size; k++) { + for(int i = 0; i < 8; i++) { + temp = els[k]*gens[i]; + + if(temp.length() > threshold /*15.*/) { + continue; + } + + is_new = words.insert(temp).second; + if(is_new == true) { + prev.push_back(temp); + } + } +} +} + +// does the axis of a given matrix go through the fundamental octagon +bool IsCanonical(const HyperbolicOctagonGroup& m) +{ + HyperbolicOctagonGroup temp = m; + + // rotate while |B1| < |B2| + SqrtField B1, B2; + SqrtField C = SqrtField(-1, -1); + for(int i = 0; i < 8 && C != C.abs(); i++) { + B1 = real(temp.B).abs(); + B2 = imag(temp.B).abs(); + C = B1 - B2; + + temp = temp.rotate(); + } + assert(C == C.abs()); + + // (2 - sqrt(2))(|B1| + (sqrt(2) - 1)|B2|) + SqrtField right = SqrtField(2, -1)*(B1 + SqrtField(-1, 1)*B2); + + // |A2| + SqrtField left = imag(temp.A).abs(); + + // left <= right -> true + C = right - left; + return C == C.abs(); +} + +void dfs(const HyperbolicOctagonGroup& m, set& visited) +{ + assert(IsCanonical(m)); + visited.insert(m); + + HyperbolicOctagonGroup candidate = m; + for(int i = 0; i < 8; i++) { + candidate = gens[i]*m*gens[(i + 4) % 8]; + if(IsCanonical(candidate) == true && visited.find(candidate) == visited.end()) { + dfs(candidate, visited); + } + } +} + +// map, m = Aux * origin * Aux^{-1} +void dfs_with_info(const pair& new_pair, + map & visited) +{ + assert(IsCanonical(new_pair.first)); + visited.insert(new_pair); + + const HyperbolicOctagonGroup& current = new_pair.first; + const HyperbolicOctagonGroup& current_factor = new_pair.second; + HyperbolicOctagonGroup candidate = current, candidate_factor = current_factor; + for(int i = 0; i < 8; i++) { + candidate = gens[i]*current*gens[(i + 4) % 8]; + if(IsCanonical(candidate) == true && visited.find(candidate) == visited.end()) { + candidate_factor = gens[i]*current_factor; + dfs_with_info(pair(candidate, candidate_factor), visited); + } + } +} + + +void dfs_with_info(const HyperbolicOctagonGroup& origin, + map& visited) +{ + assert(IsCanonical(origin)); + HyperbolicOctagonGroup id = HyperbolicOctagonGroup(Element(SqrtField(1, 0), SqrtField(0, 0)), + Element(SqrtField(0, 0), SqrtField(0, 0))); + pair new_pair(origin, id); + + dfs_with_info(new_pair, visited); +} + + +class IntersectionNumber +{ +public: + + struct Point + { + Point(double x_ = 0, double y_ = 0) : x(x_), y(y_) {} + + double x, y; + }; + + HyperbolicOctagonGroup m; + + IntersectionNumber(const HyperbolicOctagonGroup& m_) : m(m_) + {} + + long operator() () const + { + set visited; + set::iterator it, it2; + map nb_map; + map::iterator mit; + + dfs(m, visited); + + set > common; + for(it = visited.begin(); it != visited.end(); ++it) { + for(it2 = it; it2 != visited.end(); ++it2) { + if(*it == *it2) { + continue; + } + if(haveIntersection(*it, *it2) == true) { + common.clear(); + count_nb(*it, *it2, common); + + mit = nb_map.find(common.size()); + if(mit != nb_map.end()) { + mit->second += 1; + } else { + nb_map.insert(pair(common.size(), 1)); + } + } + } + } + + long nb = 0; + for(mit = nb_map.begin(); mit != nb_map.end(); mit++) { + assert( mit->second % mit->first == 0 ); + nb += mit->second/mit->first; + } + return nb; + } + + + void count_nb(const HyperbolicOctagonGroup& m1, const HyperbolicOctagonGroup& m2, set >& visited) const + { + typedef pair matrix_pair; + visited.insert(matrix_pair(m1, m2)); + + HyperbolicOctagonGroup c1 = m1, c2 = m2; + for(int i = 0; i < 8; i++) { + c1 = gens[i]*m1*gens[(i + 4) % 8]; + c2 = gens[i]*m2*gens[(i + 4) % 8]; + if(IsCanonical(c1) == true && IsCanonical(c2) == true && visited.find(matrix_pair(c1, c2)) == visited.end()) { + count_nb(c1, c2, visited); + } + } + } + +//private: + +// check whether two axis have intersection + bool haveIntersection(const HyperbolicOctagonGroup& m1, const HyperbolicOctagonGroup& m2) const + { + Point p1, p2; + intersectWithInfinity(m1, p1, p2); + + Point p3, p4; + intersectWithInfinity(m2, p3, p4); + + // orientation test + double sign1 = (p1.x - p3.x)*(p2.y - p3.y) - (p1.y - p3.y)*(p2.x - p3.x); + double sign2 = (p1.x - p4.x)*(p2.y - p4.y) - (p1.y - p4.y)*(p2.x - p4.x); + + assert( sign1 * sign2 != 0); + return (sign1 * sign2 < 0); + } + + void intersectWithInfinity(const HyperbolicOctagonGroup& m, Point& p1, Point& p2) const + { + Element a = m.A, b = m.B, factor = m.factor; + + Element four = Element(SqrtField(4, 0), SqrtField(0, 0)); + Element two = Element(SqrtField(2, 0), SqrtField(0, 0)); + + Element D = (a - conj(a))*(a - conj(a)); + D += four*b*conj(b)*factor; + Element T1 = conj(a) - a; + Element T2 = two*conj(b); + + complex d = m.toComplexDouble(D); + complex t1 = m.toComplexDouble(T1); + complex t2 = m.toComplexDouble(T2); + complex au = complex(m.factor.l + sqrt(2.)*m.factor.r, 0); + + complex z1 = (t1 + sqrt(d))/(t2*sqrt(au)); + complex z2 = (t1 - sqrt(d))/(t2*sqrt(au)); + + p1 = Point(real(z1), imag(z1)); + p2 = Point(real(z2), imag(z2)); + + assert(p1.x*p1.x + p1.y*p1.y > 0.99 && p1.x*p1.x + p1.y*p1.y < 1.01); + assert(p2.x*p2.x + p2.y*p2.y > 0.99 && p2.x*p2.x + p2.y*p2.y < 1.01); + } + + +}; + +void Delete(const set& canonical_set, vector& output) +{ + set redundant; + + set::iterator it; + for(it = canonical_set.begin(); it != canonical_set.end(); ++it) { + if(redundant.find(*it) != redundant.end()) { + continue; + } + + set visited; + dfs(*it, visited); + visited.erase(*it); + + redundant.insert(visited.begin(), visited.end()); + output.push_back(*it); + } +} + +void generate_unique_words(vector& output, double threshold = 10, int word_length = 13) +{ + get_generators(gens); + + set unique_words; + vector temp; + generate_words(unique_words, temp, word_length, threshold); + + double l = 0; + set::iterator uit; + for(uit = unique_words.begin(); uit != unique_words.end(); ++uit) { + l = uit->length(); + if(0. < l && l < threshold) { + output.push_back( *uit ); + } + } + + cout << "nb of unique words " << output.size() << endl; +} + +// words that correspond to union of 1-cycles +void generate_words_union_1_cycles(vector& out) +{ + if(gens.size() == 0) { + get_generators(gens); + } + HyperbolicOctagonGroup f[4] = {gens[0], gens[5], gens[2], gens[7]}; + + HyperbolicOctagonGroup F[8] = { + f[0]*f[0].inverse(), + f[0], + f[0]*f[1], + f[0]*f[1]*f[2], + f[0]*f[1]*f[2]*f[3], + f[3]*f[2]*f[1], + f[3]*f[2], + f[3] + }; + + long counter = 0; + for(int i = 1; i < 5; i++) { + for(int j = 0; j < 8; j++) { + for(int k = 0; k < 8; k++) { + if (j == k) { + continue; + } + // intersection + if ((0 < j && j < i) && (i < k)) { + continue; + } + if ((0 < k && k < i) && (i < j)) { + continue; + } + counter++; + + HyperbolicOctagonGroup Tr = F[i]*F[k].inverse()*F[j]; + // check that Tr != Id + if(Tr.length() > 2.) { + out.push_back(Tr); + } + } + } + } + cout << counter << endl; + +} + + + +#endif // CGAL_HYPERBOLIC_OCTAGON_GROUP_H + + + diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Octagon_matrix.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Octagon_matrix.h deleted file mode 100644 index 9dce5b1e1cf..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Octagon_matrix.h +++ /dev/null @@ -1,568 +0,0 @@ - -#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_OCTAGON_MATRIX_H -#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_OCTAGON_MATRIX_H - - -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include - - -template -class Octagon_matrix -{ -public: - typedef Octagon_matrix Self; - typedef complex Extended_field; - - Extended_field M11; - Extended_field M12; - string label; - - Octagon_matrix(const Extended_field& M11_, const Extended_field& M12_, - const string& label_ = string("") ) : - M11(M11_), M12(M12_), label(label_) {} - - string merge_labels(const Self& rh) const - { - return label + "*" + rh.label; - } - - Self operator*(const Self& rh) const - { - return Self( M11*rh.M11 + factor*M12*conj(rh.M12), - M11*rh.M12 + M12*conj(rh.M11), merge_labels(rh) ); - } - - Self inverse() const - { - Self inv = Self(conj(M11), -M12); - string inv_label; - for(long i = label.size() - 1; i >= 0; i--) { - if(label[i] >= 'a' && label[i] <= 'd') { - inv_label.push_back(label[i]); - inv_label.push_back('^'); - inv_label.push_back('-'); - inv_label.push_back('1'); - } - if(label[i] == '*') { - inv_label.push_back('*'); - } - if(label[i] == '1') { - assert(i - 3 >= 0); - assert(label[i - 3] >= 'a' && label[i - 3] <= 'd'); - inv_label.push_back(label[i - 3]); - i = i - 3; - } - } - inv.label = inv_label; - return inv; - } - - // rotation \pi/4 - Octagon_matrix rotate() const - { - Sqrt_field B1 = real(M12); - Sqrt_field B2 = imag(M12); - - // sqrt(2) - Sqrt_field k = Sqrt_field(0, 1); - - Sqrt_field BB1 = (B1 - B2)*k; - Sqrt_field BB2 = (B1 + B2)*k; - - assert(BB2.l % 2 == 0 && BB2.r % 2 == 0); - - BB1.l = BB1.l/2; - BB1.r = BB1.r/2; - BB2.l = BB2.l/2; - BB2.r = BB2.r/2; - - return Octagon_matrix(M11, Extended_field(BB1, BB2)); - } - - Sqrt_field trace() const - { - return Sqrt_field(2, 0)*real(M11); - } - - double length() const - { - typedef long double ld; - - ld l = real(M11).l; - ld r = real(M11).r; - ld tr = l + sqrt(2.)*r; - if (tr < 0) { - tr = -tr; - } - - return 2.*acosh(tr); - } - - // determinant == 1 - Extended_field det() const - { - return norm(M11) - factor * norm(M12); - } - - static complex toComplexDouble(Extended_field M) //const - { - Sqrt_field rl = real(M); - Sqrt_field img = imag(M); - - return complex(rl.l + sqrt(2.)*rl.r, img.l + sqrt(2.)*img.r); - } - - pair apply(double x, double y) - { - typedef complex Cmpl; - Cmpl m11 = toComplexDouble(M11); - Cmpl m12 = toComplexDouble(M12); - - double ax = sqrt(factor.l + sqrt(2.)*factor.r); - - Cmpl z(x, y); - Cmpl res = (m11*z + ax*m12)/(ax*(conj(m12)*z) + conj(m11)); - return pair(real(res), imag(res)); - } - -//private: - static Sqrt_field factor; -}; - - -template -Sqrt_field Octagon_matrix::factor = Sqrt_field(-1, 1); - -// just to give an order(ing) -template -bool operator < (const complex >& lh, - const complex >& rh) -{ - if (real(lh) < real(rh)) { - return true; - } - - if (real(lh) == real(rh)) { - if (imag(lh) < imag(rh)) { - return true; - } - } - - return false; -} - -// just to order octagon_matrices -template -bool operator < (const Octagon_matrix& lh, - const Octagon_matrix& rh) -{ - if (lh.M11 < rh.M11) { - return true; - } - - if (lh.M11 == rh.M11 ) { - if (lh.M12 < rh.M12) { - return true; - } - } - - return false; -} - -template -bool operator == (const Octagon_matrix& lh, - const Octagon_matrix& rh) -{ - return (lh.M11 == rh.M11 && lh.M12 == rh.M12); -} - -template -ostream& operator<<(ostream& os, const Octagon_matrix& m) -{ - os << m.M11 << " " << m.M12; - return os; -} - -typedef long long ll; -typedef Sqrt_field SqrtField; -typedef Octagon_matrix OctagonMatrix; -typedef OctagonMatrix::Extended_field Entry; - - -enum Direction { - A = 0, // 0 - InvB, // 1 - C, // 2 - InvD, // 3 - InvA, // 4 - B, // 5 - InvC, // 6 - D // 7 -}; - - -void get_generators(vector& gens) -{ - // This is a in the matrix, equal to sqrt(2) + 1 - Entry M11 = Entry(SqrtField(1, 1), SqrtField(0, 0)); - - // This vector should hold all other elements, results of the exponentials for various k - vector M12(8, Entry(SqrtField(0, 0), SqrtField(0, 0))); - - // Set everything manually - - /* - M12[A] = M11 * Entry(SqrtField(1, 0), SqrtField(-1, 0)); - M12[InvB] = M11 * Entry(SqrtField(1, 0), SqrtField(1, 0)); - M12[C] = M11 * Entry(SqrtField(-1, 0), SqrtField(1, 0)); - M12[InvD] = M11 * Entry(SqrtField(0, 0), SqrtField(0, -1)); - M12[InvA] = M11 * Entry(SqrtField(0, -1), SqrtField(0, 0)); - M12[B] = M11 * Entry(SqrtField(0, 0), SqrtField(0, 1)); - M12[InvC] = M11 * Entry(SqrtField(0, 1), SqrtField(0, 0)); - M12[D] = M11 * Entry(SqrtField(-1, 0), SqrtField(-1, 0)); - */ - - M12[C] = M11 * Entry(SqrtField(0, 1), SqrtField(0, 0)); - M12[D] = M11 * Entry(SqrtField(1, 0), SqrtField(1, 0)); - M12[InvA] = M11 * Entry(SqrtField(0, 0), SqrtField(0, 1)); - M12[InvB] = M11 * Entry(SqrtField(-1, 0), SqrtField(1, 0)); - M12[InvC] = M11 * Entry(SqrtField(0, -1), SqrtField(0, 0)); - M12[InvD] = M11 * Entry(SqrtField(-1, 0), SqrtField(-1, 0)); - M12[A] = M11 * Entry(SqrtField(0, 0), SqrtField(0, -1)); - M12[B] = M11 * Entry(SqrtField(1, 0), SqrtField(-1, 0)); - - string labels[8] = {string("a"), string("b^-1"), string("c"), string("d^-1"), - string("a^-1"), string("b"), string("c^-1"), string("d")}; - for(int i = 0; i < 8; i++) { - gens.push_back(OctagonMatrix(M11, M12[i], labels[i])); - } -} - -// a, b, c, d, a^-1, b^-1, c^-1, d^-1 -vector gens; - -bool IsCanonical(const OctagonMatrix& m); - -void generate_words( set& words, vector& prev, int depth, double threshold ) -{ - if (depth == 1) { - for(int i = 0; i < 8; i++) { - words.insert( gens[i] ); - prev.push_back( gens[i] ); - } - return; - } - - vector els; - generate_words( words, els, depth - 1, threshold); - - OctagonMatrix temp = OctagonMatrix(Entry(), Entry()); - ll size = els.size(); - bool is_new = false; - for(ll k = 0; k < size; k++) { - for(int i = 0; i < 8; i++) { - temp = els[k]*gens[i]; - - if(temp.length() > threshold /*15.*/) { - continue; - } - - is_new = words.insert(temp).second; - if(is_new == true) { - prev.push_back(temp); - } - } - } -} - -// does the axis of a given matrix go through the fundamental octagon -bool IsCanonical(const OctagonMatrix& m) -{ - OctagonMatrix temp = m; - - // rotate while |B1| < |B2| - SqrtField B1, B2; - SqrtField C = SqrtField(-1, -1); - for(int i = 0; i < 8 && C != C.abs(); i++) { - B1 = real(temp.M12).abs(); - B2 = imag(temp.M12).abs(); - C = B1 - B2; - - temp = temp.rotate(); - } - assert(C == C.abs()); - - // (2 - sqrt(2))(|B1| + (sqrt(2) - 1)|B2|) - SqrtField right = SqrtField(2, -1)*(B1 + SqrtField(-1, 1)*B2); - - // |A2| - SqrtField left = imag(temp.M11).abs(); - - // left <= right -> true - C = right - left; - return C == C.abs(); -} - -void dfs(const OctagonMatrix& m, set& visited) -{ - assert(IsCanonical(m)); - visited.insert(m); - - OctagonMatrix candidate = m; - for(int i = 0; i < 8; i++) { - candidate = gens[i]*m*gens[(i + 4) % 8]; - if(IsCanonical(candidate) == true && visited.find(candidate) == visited.end()) { - dfs(candidate, visited); - } - } -} - -// map, m = Aux * origin * Aux^{-1} -void dfs_with_info(const pair& new_pair, - map & visited) -{ - assert(IsCanonical(new_pair.first)); - visited.insert(new_pair); - - const OctagonMatrix& current = new_pair.first; - const OctagonMatrix& current_factor = new_pair.second; - OctagonMatrix candidate = current, candidate_factor = current_factor; - for(int i = 0; i < 8; i++) { - candidate = gens[i]*current*gens[(i + 4) % 8]; - if(IsCanonical(candidate) == true && visited.find(candidate) == visited.end()) { - candidate_factor = gens[i]*current_factor; - dfs_with_info(pair(candidate, candidate_factor), visited); - } - } -} - - -void dfs_with_info(const OctagonMatrix& origin, - map& visited) -{ - assert(IsCanonical(origin)); - OctagonMatrix id = OctagonMatrix(Entry(SqrtField(1, 0), SqrtField(0, 0)), - Entry(SqrtField(0, 0), SqrtField(0, 0))); - pair new_pair(origin, id); - - dfs_with_info(new_pair, visited); -} - - -class IntersectionNumber -{ -public: - -struct Point -{ - Point(double x_ = 0, double y_ = 0) : x(x_), y(y_) {} - - double x, y; -}; - -OctagonMatrix m; - -IntersectionNumber(const OctagonMatrix& m_) : m(m_) -{} - -long operator() () const -{ - set visited; - set::iterator it, it2; - map nb_map; - map::iterator mit; - - dfs(m, visited); - - set > common; - for(it = visited.begin(); it != visited.end(); ++it) { - for(it2 = it; it2 != visited.end(); ++it2) { - if(*it == *it2) { - continue; - } - if(haveIntersection(*it, *it2) == true) { - common.clear(); - count_nb(*it, *it2, common); - - mit = nb_map.find(common.size()); - if(mit != nb_map.end()) { - mit->second += 1; - } else { - nb_map.insert(pair(common.size(), 1)); - } - } - } - } - - long nb = 0; - for(mit = nb_map.begin(); mit != nb_map.end(); mit++) { - assert( mit->second % mit->first == 0 ); - nb += mit->second/mit->first; - } - return nb; -} - - -void count_nb(const OctagonMatrix& m1, const OctagonMatrix& m2, set >& visited) const -{ - typedef pair matrix_pair; - visited.insert(matrix_pair(m1, m2)); - - OctagonMatrix c1 = m1, c2 = m2; - for(int i = 0; i < 8; i++) { - c1 = gens[i]*m1*gens[(i + 4) % 8]; - c2 = gens[i]*m2*gens[(i + 4) % 8]; - if(IsCanonical(c1) == true && IsCanonical(c2) == true && visited.find(matrix_pair(c1, c2)) == visited.end()) { - count_nb(c1, c2, visited); - } - } -} - -//private: - -// check whether two axis have intersection -bool haveIntersection(const OctagonMatrix& m1, const OctagonMatrix& m2) const -{ - Point p1, p2; - intersectWithInfinity(m1, p1, p2); - - Point p3, p4; - intersectWithInfinity(m2, p3, p4); - - // orientation test - double sign1 = (p1.x - p3.x)*(p2.y - p3.y) - (p1.y - p3.y)*(p2.x - p3.x); - double sign2 = (p1.x - p4.x)*(p2.y - p4.y) - (p1.y - p4.y)*(p2.x - p4.x); - - assert( sign1 * sign2 != 0); - return (sign1 * sign2 < 0); -} - -void intersectWithInfinity(const OctagonMatrix& m, Point& p1, Point& p2) const -{ - Entry a = m.M11, b = m.M12, factor = m.factor; - - Entry four = Entry(SqrtField(4, 0), SqrtField(0, 0)); - Entry two = Entry(SqrtField(2, 0), SqrtField(0, 0)); - - Entry D = (a - conj(a))*(a - conj(a)); - D += four*b*conj(b)*factor; - Entry T1 = conj(a) - a; - Entry T2 = two*conj(b); - - complex d = m.toComplexDouble(D); - complex t1 = m.toComplexDouble(T1); - complex t2 = m.toComplexDouble(T2); - complex au = complex(m.factor.l + sqrt(2.)*m.factor.r, 0); - - complex z1 = (t1 + sqrt(d))/(t2*sqrt(au)); - complex z2 = (t1 - sqrt(d))/(t2*sqrt(au)); - - p1 = Point(real(z1), imag(z1)); - p2 = Point(real(z2), imag(z2)); - - assert(p1.x*p1.x + p1.y*p1.y > 0.99 && p1.x*p1.x + p1.y*p1.y < 1.01); - assert(p2.x*p2.x + p2.y*p2.y > 0.99 && p2.x*p2.x + p2.y*p2.y < 1.01); -} - - -}; - -void Delete(const set& canonical_set, vector& output) -{ - set redundant; - - set::iterator it; - for(it = canonical_set.begin(); it != canonical_set.end(); ++it) { - if(redundant.find(*it) != redundant.end()) { - continue; - } - - set visited; - dfs(*it, visited); - visited.erase(*it); - - redundant.insert(visited.begin(), visited.end()); - output.push_back(*it); - } -} - -void generate_unique_words(vector& output, double threshold = 10, int word_length = 13) -{ - get_generators(gens); - - set unique_words; - vector temp; - generate_words(unique_words, temp, word_length, threshold); - - double l = 0; - set::iterator uit; - for(uit = unique_words.begin(); uit != unique_words.end(); ++uit) { - l = uit->length(); - if(0. < l && l < threshold) { - output.push_back( *uit ); - } - } - - cout << "nb of unique words " << output.size() << endl; -} - -// words that correspond to union of 1-cycles -void generate_words_union_1_cycles(vector& out) -{ - if(gens.size() == 0) { - get_generators(gens); - } - OctagonMatrix f[4] = {gens[0], gens[5], gens[2], gens[7]}; - - OctagonMatrix F[8] = { - f[0]*f[0].inverse(), - f[0], - f[0]*f[1], - f[0]*f[1]*f[2], - f[0]*f[1]*f[2]*f[3], - f[3]*f[2]*f[1], - f[3]*f[2], - f[3] - }; - - long counter = 0; - for(int i = 1; i < 5; i++) { - for(int j = 0; j < 8; j++) { - for(int k = 0; k < 8; k++) { - if (j == k) { - continue; - } - // intersection - if ((0 < j && j < i) && (i < k)) { - continue; - } - if ((0 < k && k < i) && (i < j)) { - continue; - } - counter++; - - OctagonMatrix Tr = F[i]*F[k].inverse()*F[j]; - // check that Tr != Id - if(Tr.length() > 2.) { - out.push_back(Tr); - } - } - } - } - cout << counter << endl; - -} - -#endif - - - diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h index efa9186df32..4b83790eecd 100644 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h @@ -22,9 +22,9 @@ #include -// I needed to use a non-reserved word for this thing, there were problems when I (stupidly) used TYPE +// I needed to use a non-reserved word for this thing, there were problems when I (iiordano) (stupidly) used TYPE #define CLEARLY_MY_TYPE 2 // 1 = Cross, 2 = Diametric -#define USE_TEST_THINGS 1 // 1 = True, 0 = False +#define USE_TEST_THINGS 0 // 1 = True, 0 = False #if USE_TEST_THINGS == 1 #include diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_construct_hyperbolic_point_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_construct_hyperbolic_point_2.h new file mode 100644 index 00000000000..18cf0c36200 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_construct_hyperbolic_point_2.h @@ -0,0 +1,66 @@ +// Copyright (c) 1999-2004,2006-2009,2014-2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Monique Teillaud +// Sylvain Pion +// Andreas Fabri +// Nico Kruithof +// Manuel Caroli +// Aymeric Pellé +#ifndef CGAL_PERIODIC_4_CONSTRUCT_HYPERBOLIC_POINT_2_H +#define CGAL_PERIODIC_4_CONSTRUCT_HYPERBOLIC_POINT_2_H + +#include +#include + +namespace CGAL +{ +template < typename K, typename Construct_point_2_base> +class Periodic_4_construct_hyperbolic_point_2 : public Construct_point_2_base +{ + typedef K Kernel; + +public: + typedef typename Kernel::Point_2 Point; + typedef HyperbolicOctagonGroup Offset; + typedef typename Kernel::Circle_2 Circle_2; + + typedef Point result_type; + + Periodic_4_construct_hyperbolic_point_2(const Circle_2 & dom) : _dom(dom) { } + + Point operator() ( const Point& p, const Offset& o ) const { + /* + return Point( p.x()+(_dom.xmax()-_dom.xmin())*o.x(), + p.y()+(_dom.ymax()-_dom.ymin())*o.y(), + p.z()+(_dom.zmax()-_dom.zmin())*o.z()); */ + pair pt = o.apply(to_double(p.x()), to_double(p.y())); + Point r = Point(pt.first, pt.second); + + // Do we need this? I (iiordano) think it's a good idea... + assert( _dom.has_on_bounded_side(r) ); + + return r; + } + +private: + Circle_2 _dom; +}; +} + +#endif diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h.misha b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h.misha deleted file mode 100644 index ad5af73f686..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h.misha +++ /dev/null @@ -1,570 +0,0 @@ -// Copyright (c) 2010 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org); you may redistribute it under -// the terms of the Q Public License version 1.0. -// See the file LICENSE.QPL distributed with CGAL. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/candidate-packages/Triangulation_2/include/CGAL/Delaunay_triangulation_2.h $ -// $Id: Delaunay_triangulation_2.h 57509 2010-07-15 09:14:09Z sloriot $ -// -// -// Author(s) : Mikhail Bogdanov - -#ifndef CGAL_DELAUNAY_HYPERBOLIC_TRIANGULATION_2_H -#define CGAL_DELAUNAY_HYPERBOLIC_TRIANGULATION_2_H - -#include -#include -#include - -#include -#include - -namespace CGAL { - -template < class Gt, - class Tds = Triangulation_data_structure_2 < - Triangulation_vertex_base_2, - Triangulation_face_base_with_info_2 > > -class Delaunay_hyperbolic_triangulation_2 : public Delaunay_triangulation_2 -{ -public: - typedef Delaunay_hyperbolic_triangulation_2 Self; - typedef Delaunay_triangulation_2 Base; - - typedef Triangulation_face_base_with_info_2 Face_base; - typedef typename Face_base::Info Face_info; - typedef typename Base::size_type size_type; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Face_handle Face_handle; - typedef typename Base::Edge Edge; - -#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 - using Base::cw; - using Base::ccw; - using Base::geom_traits; -#endif - - typedef typename Base::Edge_circulator Edge_circulator; - typedef typename Base::Face_circulator Face_circulator; - typedef typename Base::Vertex_circulator Vertex_circulator; - - typedef typename Base::All_vertices_iterator All_vertices_iterator; - typedef typename Base::All_edges_iterator All_edges_iterator; - typedef typename Base::All_faces_iterator All_faces_iterator; - - typedef Gt Geom_traits; - typedef typename Geom_traits::FT FT; - typedef typename Geom_traits::Point_2 Point; - typedef typename Geom_traits::Segment_2 Segment; - - /* -#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 - using Triangulation::side_of_oriented_circle; - using Triangulation::circumcenter; - using Triangulation::collinear_between; - using Triangulation::test_dim_down; - using Triangulation::make_hole; - using Triangulation::fill_hole_delaunay; - using Triangulation::delete_vertex; -#endif -*/ - - Delaunay_hyperbolic_triangulation_2(const Gt& gt = Gt()) - : Delaunay_triangulation_2(gt) {} - - Delaunay_hyperbolic_triangulation_2( - const Delaunay_hyperbolic_triangulation_2 &tr) - : Delaunay_triangulation_2(tr) - { CGAL_triangulation_postcondition( this->is_valid() );} - - - void mark_star(Vertex_handle v) const - { - if(!is_star_bounded(v)) { - mark_star_faces(v); - } - } - - Vertex_handle insert(const Point &p, - Face_handle start = Face_handle() ) - { - Vertex_handle v = Base::insert(p, start); - mark_star(v); - - return v; - } - - Vertex_handle insert(const Point& p, - typename Base::Locate_type lt, - Face_handle loc, int li ) - { - Vertex_handle v = Base::insert(p, lt, loc, li); - mark_star(v); - - return v; - } - -#ifndef CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - template < class InputIterator > - std::ptrdiff_t - insert( InputIterator first, InputIterator last, - typename boost::enable_if< - boost::is_base_of< - Point, - typename std::iterator_traits::value_type - > - >::type* = NULL - ) -#else - template < class InputIterator > - std::ptrdiff_t - insert(InputIterator first, InputIterator last) -#endif //CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - { - size_type n = Base::insert(first, last); - - mark_faces(); - - return n; - } - - //test version of insert function - -#ifndef CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - template < class InputIterator > - std::ptrdiff_t - insert2( InputIterator first, InputIterator last, - typename boost::enable_if< - boost::is_base_of< - Point, - typename std::iterator_traits::value_type - > - >::type* = NULL - ) -#else - template < class InputIterator > - std::ptrdiff_t - insert_2(InputIterator first, InputIterator last) -#endif //CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO - { - size_type n = this->number_of_vertices(); - - spatial_sort(first, last, geom_traits()); - Face_handle f; - while(first != last) { - f = insert (*first++, f)->face(); - } - - return this->number_of_vertices() - n; - } - - bool is_infinite(Vertex_handle v) const - { - return Base::is_infinite(v); - } - - bool is_infinite(Face_handle f) const - { - return has_infinite_vertex(f) || is_finite_invisible(f); - } - - bool is_infinite(Face_handle f, int i) const - { - return has_infinite_vertex(f, i) || is_finite_invisible(f, i); - } - - bool is_infinite(const Edge& e) const - { - return is_infinite(e.first, e.second); - } - - bool is_infinite(const Edge_circulator& ec) const - { - return is_infinite(*ec); - } - - bool is_infinite(const All_edges_iterator& ei) const - { - return is_infinite(*ei); - } - -private: - - bool has_infinite_vertex(Face_handle f) const - { - return Base::is_infinite(f); - } - - bool has_infinite_vertex(Face_handle f, int i) const - { - return Base::is_infinite(f, i); - } - - bool has_infinite_vertex(const Edge& e) const - { - return Base::is_infinite(e); - } - - int get_finite_invisible_edge(Face_handle f) const - { - assert(is_finite_invisible(f)); - - return f->info().get_invisible_edge(); - } - - bool is_finite_invisible(Face_handle f) const - { - return f->info().is_finite_invisible(); - } - - bool is_finite_invisible(Face_handle f, int i) const - { - if(this->dimension() <= 1) { - return false; - } - - if(is_finite_invisible(f) && get_finite_invisible_edge(f) == i) { - return true; - } - - // another incident face and corresponding index - Face_handle f2 = f->neighbor(i); - int i2 = f2->index(f); - - if(is_finite_invisible(f2) && get_finite_invisible_edge(f2) == i2) { - return true; - } - - return false; - } - - bool is_finite_invisible(const Edge& e) const - { - return is_finite_invisible(e.first, e.second); - } - - // Depth-first search (dfs) and marking the finite invisible faces. - void mark_faces() const - { - if(this->dimension() <= 1) return; - - std::set visited_faces; - - // maintain a stack to be able to backtrack - // to the most recent faces which neighbors are not visited - std::stack backtrack; - - // start from a face with infinite vertex - Face_handle current = Base::infinite_face(); - - // mark it as visited - visited_faces.insert(current); - - // put the element whose neighbors we are going to explore. - backtrack.push(current); - - // test whether a face is finite invisible or not - Mark_face test(*this); - - Face_handle next; - Face_info face_info; - - while(!backtrack.empty()) { - // take a face - current = backtrack.top(); - - // start visiting the neighbors - int i = 0; - for(; i < 3; i++) { - next = current->neighbor(i); - - // if a neighbor is already visited, then stop going deeper - if(visited_faces.find(next) != visited_faces.end()) { - continue; - } - - visited_faces.insert(next); - mark_face(next, test); - - // go deeper if the neighbor is infinite - if(is_infinite(next)) { - backtrack.push(next); - break; - } - } - - // if all the neighbors are already visited, then remove "current" face. - if(i == 3) { - backtrack.pop(); - } - } - - } - - // check if a star is bounded by finite faces - // TODO: rename this function name - bool is_star_bounded(Vertex_handle v) const - { - if(this->dimension() <= 1) { - return true; - } - - Face_handle f = v->face(); - Face_handle next; - int i; - Face_handle start(f); - - Face_handle opposite_face; - - do { - i = f->index(v); - next = f->neighbor(ccw(i)); // turn ccw around v - - opposite_face = f->neighbor(i); - if(this->is_infinite(opposite_face)) { - return false; - } - - f = next; - } while(next != start); - - return true; - } - - //use the function: insert_and_give_new_faces? - - void mark_star_faces(Vertex_handle v) const - { - // TODO: think of it - if(this->dimension() <= 1) return; - - Mark_face test(*this); - - Face_handle f = v->face(); - Face_handle start(f), next; - int i; - do { - i = f->index(v); - next = f->neighbor(ccw(i)); // turn ccw around v - - mark_face(f, test); - - f = next; - } while(next != start); - return; - } - - template - void mark_face(const Face_handle& f, const Mark_face_test& test) const - { - f->info() = test(f); - } - - void mark_face(const Face_handle& f) const - { - Mark_face test(*this); - mark_face(f, test); - } - - class Mark_face - { - public: - Mark_face(const Self& tr) : - _tr(tr) - {} - - Face_info operator ()(const Face_handle& f) const - { - typedef typename Gt::Is_hyperbolic Is_hyperbolic; - - Face_info info; - if(_tr.has_infinite_vertex(f)) { - return info; - } - - Point p0 = f->vertex(0)->point(); - Point p1 = f->vertex(1)->point(); - Point p2 = f->vertex(2)->point(); - int ind = 0; - - Is_hyperbolic is_hyperbolic = _tr.geom_traits().Is_hyperbolic_object(); - if(is_hyperbolic(p0, p1, p2, ind) == false) { - - info.set_finite_invisible(true); - info.set_invisible_edge(ind); - - return info; - } - - return info; - } - - private: - - Mark_face(const Mark_face&); - Mark_face& operator= (const Mark_face&); - - const Self& _tr; - }; - -public: - // This class is used to generate the Finite_*_iterators. - class Infinite_hyperbolic_tester - { - const Self *t; - public: - Infinite_hyperbolic_tester() {} - Infinite_hyperbolic_tester(const Self *tr) : t(tr) {} - - bool operator()(const All_vertices_iterator & vit) const { - return t->is_infinite(vit); - } - bool operator()(const All_faces_iterator & fit) const { - return t->is_infinite(fit); - } - bool operator()(const All_edges_iterator & eit ) const { - return t->is_infinite(eit); - } - }; - - Infinite_hyperbolic_tester - infinite_hyperbolic_tester() const - { - return Infinite_hyperbolic_tester(this); - } - - //Finite faces iterator - - class Finite_faces_iterator - : public Filter_iterator - { - typedef Filter_iterator Base; - typedef Finite_faces_iterator Self; - public: - Finite_faces_iterator() : Base() {} - Finite_faces_iterator(const Base &b) : Base(b) {} - Self & operator++() { Base::operator++(); return *this; } - Self & operator--() { Base::operator--(); return *this; } - Self operator++(int) { Self tmp(*this); ++(*this); return tmp; } - Self operator--(int) { Self tmp(*this); --(*this); return tmp; } - operator const Face_handle() const { return Base::base(); } - }; - - Finite_faces_iterator - finite_faces_begin() const - { - if ( this->dimension() < 2 ) - return finite_faces_end(); - return CGAL::filter_iterator(this->all_faces_end(), - Infinite_hyperbolic_tester(this), - this->all_faces_begin() ); - } - - Finite_faces_iterator - finite_faces_end() const - { - return CGAL::filter_iterator(this->all_faces_end(), - Infinite_hyperbolic_tester(this) ); - } - - //Finite edges iterator - - typedef Filter_iterator Finite_edges_iterator; - - Finite_edges_iterator - finite_edges_begin() const - { - if ( this->dimension() < 1 ) - return finite_edges_end(); - return CGAL::filter_iterator(this->all_edges_end(), - infinite_hyperbolic_tester(), - this->all_edges_begin()); - } - - Finite_edges_iterator - finite_edges_end() const - { - return CGAL::filter_iterator(this->all_edges_end(), - infinite_hyperbolic_tester() ); - } - - using Base::dual; - - Object - dual(const Finite_edges_iterator& ei) const - { - return this->dual(*ei); - } - - Object - dual(const Edge &e) const - { - CGAL_triangulation_precondition (!this->is_infinite(e)); - - if(this->dimension() == 1) { - const Point& p = (e.first)->vertex(cw(e.second))->point(); - const Point& q = (e.first)->vertex(ccw(e.second))->point(); - - // hyperbolic line - Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - return make_object(line); - } - - // incident faces - Face_handle f1 = e.first; - int i1 = e.second; - - Face_handle f2 = f1->neighbor(i1); - int i2 = f2->index(f1); - - // boths faces are infinite, but the incident edge is finite - if(is_infinite(f1) && is_infinite(f2)){ - const Point& p = (f1)->vertex(cw(i1))->point(); - const Point& q = (f1)->vertex(ccw(i1))->point(); - - // hyperbolic line - Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - return make_object(line); - } - - // both faces are finite - if(!is_infinite(f1) && !is_infinite(f2)) { - - Segment s = this->geom_traits().construct_segment_2_object()(dual(f1),dual(f2)); - - return make_object(s); - } - - // one of the incident faces is infinite - Face_handle finite_face = f1; - int i = i1; - - if(is_infinite(f1)) { - finite_face = f2; - i = i2; - } - - const Point& p = finite_face->vertex(cw(i))->point(); - const Point& q = finite_face->vertex(ccw(i))->point(); - - // ToDo: Line or Segment? - // hyperbolic line and ray - Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q); - Segment ray = this->geom_traits().construct_ray_2_object()(dual(finite_face), line); - return make_object(ray); - } -}; - -} //namespace CGAL - -#endif // CGAL_DELAUNAY_HYPERBOLIC_TRIANGULATION_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h new file mode 100644 index 00000000000..5cf187620c2 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h @@ -0,0 +1,96 @@ +// Copyright (c) 2006-2009 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Nico Kruithof +// Manuel Caroli + +#ifndef CGAL_PERIODIC_4_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H +#define CGAL_PERIODIC_4_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H + +//#include +#include +#include +#include +#include + +namespace CGAL { + +template < class THT2 > +class Periodic_4_hyperbolic_Delaunay_triangulation_traits_base_2 + : public Periodic_3_triangulation_traits_base_3 +{ +public: + typedef THT2 K; + typedef HyperbolicOctagonGroup Offset; + typedef Periodic_4_hyperbolic_triangulation_traits_base_2< K > Base; + typedef Periodic_4_hyperbolic_Delaunay_triangulation_traits_base_2< K > Self; + + typedef typename Base::RT RT; + typedef typename Base::FT FT; + typedef typename Base::Point_2 Point_2; + typedef typename Base::Vector_2 Vector_2; + //typedef typename Base::Periodic_4_offset_2 Periodic_4_offset_2; + + // The next typedef is there for backward compatibility + // Some users take their point type from the traits class. + // Before this type was Point + typedef Point_2 Point; + + typedef typename Base::Segment_2 Segment_2; + typedef typename Base::Triangle_2 Triangle_2; + + // Delaunay specific predicates + typedef Traits_with_offsets_adaptor + Side_of_oriented_circle_2; + typedef Traits_with_offsets_adaptor + Compare_distance_2; + typedef Traits_with_offsets_adaptor + Side_of_bounded_disc_2; + + // Delaunay specific constructions + typedef Traits_with_hyperbolic_offsets_adaptor + Construct_circumcenter_2; + + // Operations + Side_of_oriented_circle_2 + side_of_oriented_circle_2_object() const { + return Side_of_oriented_circle_2(&this->_domain); + } + Compare_distance_2 + compare_distance_2_object() const { + return Compare_distance_2(&this->_domain); + } + Side_of_bounded_disc_2 + side_of_bounded_disc_2_object() const { + return Side_of_bounded_disc_2(&this->_domain); + } + Construct_circumcenter_2 + construct_circumcenter_2_object() const { + return Construct_circumcenter_2(&this->_domain); + } +}; + +template < typename K > +class Periodic_4_hyperbolic_Delaunay_triangulation_traits_2; + +} //namespace CGAL + +#endif // CGAL_PERIODIC_4_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_offset_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_offset_2.h deleted file mode 100644 index 5d18e629c49..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_offset_2.h +++ /dev/null @@ -1,182 +0,0 @@ -// Copyright (c) 1997-2013 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Nico Kruithof - -#ifndef CGAL_PERIODIC_4_HYPERBOLIC_OFFSET_2_H -#define CGAL_PERIODIC_4_HYPERBOLIC_OFFSET_2_H - -#include -#include -#include - -namespace CGAL -{ - -/// The class Periodic_4_hyperbolic_offset_2 is a model of the concept Periodic_2Offset_2. -class Periodic_4_hyperbolic_offset_2 -{ - // template - // friend std::ostream & operator<<(std::ostream &os, - // const Periodic_4_hyperbolic_offset_2 &off); - -public: - /// Default constructor. - Periodic_4_hyperbolic_offset_2() : _offx(0), _offy(0) {} - /// Constructs the offset (x,y). - Periodic_4_hyperbolic_offset_2(int x, int y) : _offx(x), _offy(y) {} - - /// Returns true if o is equal to (0,0). - inline bool is_null() const - { - return is_zero(); - } - /// Returns true if o is equal to (0,0). - inline bool is_zero() const - { - return ((_offx | _offy) == 0); - } - - /// Return the x-entry of o. - int& x() - { - return _offx; - } - /// Return the x-entry of o. - int x() const - { - return _offx; - } - /// Return the y-entry of o. - int& y() - { - return _offy; - } - /// Return the y-entry of o. - int y() const - { - return _offy; - } - - /// Return the i-th entry of o. - int &operator[](int i) - { - if (i == 0) return _offx; - CGAL_triangulation_assertion(i == 1); - return _offy; - } - /// Return the i-th entry of o. - int operator[](int i) const - { - if (i == 0) return _offx; - CGAL_triangulation_assertion(i == 1); - return _offy; - } - /// Add o' to o using vector addition. - void operator+=(const Periodic_4_hyperbolic_offset_2 &other) - { - _offx += other._offx; - _offy += other._offy; - } - /// Subtract o' from o using vector subtraction. - void operator-=(const Periodic_4_hyperbolic_offset_2 &other) - { - _offx -= other._offx; - _offy -= other._offy; - } - /// Return the negative vector of o. - Periodic_4_hyperbolic_offset_2 operator-() const - { - return Periodic_4_hyperbolic_offset_2(-_offx, -_offy); - } - /// Return true if o' and o represent the same vector. - bool operator==(const Periodic_4_hyperbolic_offset_2 &other) const - { - return ((_offx == other._offx) && - (_offy == other._offy)); - } - /// Return true if o' and o do not represent the same vector. - bool operator!=(const Periodic_4_hyperbolic_offset_2 &other) const - { - return ((_offx != other._offx) || - (_offy != other._offy)); - } - /// Compare o and o' lexicographically. - bool operator<(const Periodic_4_hyperbolic_offset_2 &other) const - { - if (_offx != other._offx) - return (_offx < other._offx); - else - { - return (_offy < other._offy); - } - } - - /// Return the vector sum of o and o'. - Periodic_4_hyperbolic_offset_2 operator+(const Periodic_4_hyperbolic_offset_2 &off2) const - { - return Periodic_4_hyperbolic_offset_2(_offx + off2.x(), _offy + off2.y()); - } - /// Return the vector difference of o and o'. - Periodic_4_hyperbolic_offset_2 operator-(const Periodic_4_hyperbolic_offset_2 &off2) const - { - return Periodic_4_hyperbolic_offset_2(_offx - off2.x(), _offy - off2.y()); - } - -private: - int _offx, _offy; -}; - -template -inline typename K::Point_2 operator+(const typename K::Point_2 &p, const Periodic_4_hyperbolic_offset_2 &off) -{ - return (off.is_null() ? p : typename K::Point_2(p.x() + off.x(), p.y() + off.y())); -} - -/// Inputs an Periodic_4_hyperbolic_offset_2 from is. -inline std::ostream -&operator<<(std::ostream &os, const Periodic_4_hyperbolic_offset_2 &off) -{ - if (is_ascii(os)) - os << off.x() << " " << off.y(); - else - { - write(os, off.x()); - write(os, off.y()); - } - return os; -} - -/// Outputs an Periodic_4_hyperbolic_offset_2 to os. -inline std::istream -&operator>>(std::istream &is, Periodic_4_hyperbolic_offset_2 &off) -{ - int x = 0, y = 0; - if (is_ascii(is)) - is >> x >> y; - else - { - read(is, x); - read(is, y); - } - off = Periodic_4_hyperbolic_offset_2(x, y); - return is; -} - -} //namespace CGAL - -#endif // CGAL_PERIODIC_4_HYPERBOLIC_OFFSET_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h new file mode 100644 index 00000000000..ef3396dd86f --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h @@ -0,0 +1,123 @@ +// Copyright (c) 2006-2009 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Nico Kruithof +// Manuel Caroli + +#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H +#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H + + +#include +#include +#include +#include +#include + + +namespace CGAL { + +template < class THT2 > +class Periodic_4_hyperbolic_triangulation_traits_base_2 + : public THT2 +{ +public: + typedef THT2 K; + typedef HyperbolicOctagonGroup Offset; + typedef Periodic_4_hyperbolic_triangulation_traits_base_2< K > Self; + + typedef typename K::RT RT; + typedef typename K::FT FT; + typedef typename K::Point_2 Point_2; + //typedef typename K::Vector_3 Vector_3; + typedef Offset Periodic_4_hyperbolic_offset_2; + typedef typename K::Circle_2 Circle_2; + + // The next typedef is there for backward compatibility + // Some users take their point type from the traits class. + // Before this type was Point + typedef Point_2 Point; + + + // Triangulation predicates + typedef Traits_with_hyperbolic_offsets_adaptor + Compare_xy_2; + typedef Traits_with_hyperbolic_offsets_adaptor + Orientation_2; + + // Triangulation constructions + typedef Periodic_4_construct_hyperbolic_point_2 + Construct_hyperbolic_point_2; + typedef Traits_with_hyperbolic_offsets_adaptor + Construct_segment_2; + typedef Traits_with_hyperbolic_offsets_adaptor + Construct_hyperbolic_triangle_2; + + // Access + void set_domain(const Circle_2& domain) { + _domain = domain; + } + + Circle_2 get_domain() const { + return _domain; + } + + // Operations + Compare_xy_2 + compare_xy_2_object() const { + return Compare_xy_2(&_domain); + } + Orientation_2 + orientation_2_object() const { + return Orientation_2(&_domain); + } + + Construct_hyperbolic_point_2 + construct_hyperbolic_point_2_object() const { + return Construct_point_2(_domain); + } + Construct_segment_2 + construct_segment_2_object() const { + return Construct_segment_2(&_domain); + } + + +protected: + Circle_2 _domain; +}; + +namespace future_release +{ +template < typename K > +class Periodic_4_hyperbolic_triangulation_traits_2; +} + +} //namespace CGAL + + +#include + + +namespace CGAL +{ +template < typename K > +class Periodic_4_hyperbolic_Delaunay_triangulation_traits_2; + +} + +#endif // CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Sqrt_field.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Sqrt_field.h deleted file mode 100644 index 6c7ef754025..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Sqrt_field.h +++ /dev/null @@ -1,217 +0,0 @@ -// to compile with boost headers -// g++ main.cpp -I /opt/local/include -o main - -#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQRT_FIELD_H -#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQRT_FIELD_H - -#include -#include -#include -#include -#include -#include -#include -//#include -#include -#include - -//#include - -using namespace std; - -template -class Sqrt_field -{ -public: - typedef Int Integer; - - Int l; - Int r; - -public: - Sqrt_field(Int l_ = 0, Int r_ = 0) : l(l_), r(r_) {} - - Sqrt_field operator + (const Sqrt_field& rh) const - { - Sqrt_field temp = *this; - temp += rh; - - return temp; - } - - Sqrt_field& operator += (const Sqrt_field& rh) - { - l += rh.l; - r += rh.r; - return *this; - } - - Sqrt_field& operator -= (const Sqrt_field& rh) - { - l -= rh.l; - r -= rh.r; - return *this; - } - - Sqrt_field operator - (const Sqrt_field& rh) const - { - return Sqrt_field( l - rh.l, r - rh.r ); - } - - Sqrt_field operator - () const - { - return Sqrt_field(-l, -r); - } - - Sqrt_field operator * (const Sqrt_field& rh) const - { - Sqrt_field temp = *this; - temp *= rh; - - return temp; - } - - Sqrt_field& operator *= (const Sqrt_field& rh) - { - Int l1 = l * rh.l + 2*r*rh.r; - Int r1 = r*rh.l + l*rh.r; - l = l1; - r = r1; - - return *this; - } - - Sqrt_field abs() const - { - if (l >= 0 && r >= 0) { - return *this; - } - - if (l <= 0 && r <= 0) { - return Sqrt_field(-l, -r); - } - - if(l*l*l - 2*r*r*l >= 0) { - return *this; - } - return Sqrt_field(-l, -r); - } -}; - -template -bool operator < (const Sqrt_field& lh, const Sqrt_field& rh) -{ - if( (lh.l - rh.l)*(lh.l - rh.l) > 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((lh.l - rh.l) < 0) { - return true; - } - } - - if( (lh.l - rh.l)*(lh.l - rh.l) < 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((rh.r - lh.r) >= 0) { - return true; - } - } - return false; - - /* - if (lh.l < rh.l) { - return true; - } - - if (lh.l == rh.l) { - if (lh.r < rh.r) { - return true; - } - }*/ - - return false; -} -/*Seems correct! -template -bool operator > (const Sqrt_field& lh, const Sqrt_field& rh) -{ - if( (lh.l - rh.l)*(lh.l - rh.l) > 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((lh.l - rh.l) >= 0 && (rh.r - lh.r) >= 0) { - return true; - } - } - - if( (lh.l - rh.l)*(lh.l - rh.l) < 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((lh.l - rh.l) <= 0 && (rh.r - lh.r) <= 0) { - return true; - } - } - return false; -}*/ - -template -bool operator >= (const Sqrt_field& lh, const Sqrt_field& rh) -{ - if( (lh.l - rh.l)*(lh.l - rh.l) >= 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((lh.l - rh.l) >= 0) { - return true; - } - } - - if( (lh.l - rh.l)*(lh.l - rh.l) < 2*(rh.r - lh.r)*(rh.r - lh.r) ) { - if ((rh.r - lh.r) < 0) { - return true; - } - } - return false; -} - -template -bool operator == (const Sqrt_field& lh, const Sqrt_field& rh) -{ - return (lh.l == rh.l && lh.r == rh.r); -} - -template -bool operator != (const Sqrt_field& lh, const Sqrt_field& rh) -{ - return !(lh == rh); -} - -template -Sqrt_field operator * (const Int& val, const Sqrt_field& rh) -{ - Sqrt_field temp = rh; - temp.l *= val; - temp.r *= val; - - return temp; -} - -template -ostream& operator<<(ostream& os, const Sqrt_field& nb) -{ - os << nb.l << " " << nb.r; - return os; -} - -// begin MT - added to please libc++ -template -bool isnan(const Sqrt_field& x) -{ - return false; // MT TBD -} - -template -bool isinf(const Sqrt_field& x) -{ - return false; // MT TBD -} - -template -Sqrt_field copysign(const Sqrt_field& x, const Sqrt_field& y) -{ - if (y.abs()==y) - { return x;} - return -x; -} -// end MT - added to please libc++ - -#endif // CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQRT_FIELD_H - diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Square_root_2_field.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Square_root_2_field.h new file mode 100644 index 00000000000..204042146c8 --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Square_root_2_field.h @@ -0,0 +1,199 @@ +// Copyright (c) 2010 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// +// Author(s) : Mikhail Bogdanov + +#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQUARE_ROOT_2_FIELD_H +#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQUARE_ROOT_2_FIELD_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +template +class Square_root_2_field +{ + +private: + typedef Int Integer; + +public: + Int l; + Int r; + + Square_root_2_field(Int l_ = 0, Int r_ = 0) : l(l_), r(r_) {} + + Square_root_2_field operator + (const Square_root_2_field& rh) const + { + Square_root_2_field temp = *this; + temp += rh; + + return temp; + } + + Square_root_2_field& operator += (const Square_root_2_field& rh) + { + l += rh.l; + r += rh.r; + return *this; + } + + Square_root_2_field& operator -= (const Square_root_2_field& rh) + { + l -= rh.l; + r -= rh.r; + return *this; + } + + Square_root_2_field operator - (const Square_root_2_field& rh) const + { + return Square_root_2_field( l - rh.l, r - rh.r ); + } + + Square_root_2_field operator - () const + { + return Square_root_2_field(-l, -r); + } + + Square_root_2_field operator * (const Square_root_2_field& rh) const + { + Square_root_2_field temp = *this; + temp *= rh; + + return temp; + } + + Square_root_2_field& operator *= (const Square_root_2_field& rh) + { + Int l1 = l * rh.l + 2*r*rh.r; + Int r1 = r*rh.l + l*rh.r; + l = l1; + r = r1; + + return *this; + } + + Square_root_2_field abs() const + { + if (l >= 0 && r >= 0) { + return *this; + } + + if (l <= 0 && r <= 0) { + return Square_root_2_field(-l, -r); + } + + if(l*l*l - 2*r*r*l >= 0) { + return *this; + } + return Square_root_2_field(-l, -r); + } +}; + +template +bool operator < (const Square_root_2_field& lh, const Square_root_2_field& rh) +{ + if( (lh.l - rh.l)*(lh.l - rh.l) > 2*(rh.r - lh.r)*(rh.r - lh.r) ) { + if ((lh.l - rh.l) < 0) { + return true; + } + } + + if( (lh.l - rh.l)*(lh.l - rh.l) < 2*(rh.r - lh.r)*(rh.r - lh.r) ) { + if ((rh.r - lh.r) >= 0) { + return true; + } + } + return false; +} + + +template +bool operator >= (const Square_root_2_field& lh, const Square_root_2_field& rh) +{ + if( (lh.l - rh.l)*(lh.l - rh.l) >= 2*(rh.r - lh.r)*(rh.r - lh.r) ) { + if ((lh.l - rh.l) >= 0) { + return true; + } + } + + if( (lh.l - rh.l)*(lh.l - rh.l) < 2*(rh.r - lh.r)*(rh.r - lh.r) ) { + if ((rh.r - lh.r) < 0) { + return true; + } + } + return false; +} + +template +bool operator == (const Square_root_2_field& lh, const Square_root_2_field& rh) +{ + return (lh.l == rh.l && lh.r == rh.r); +} + +template +bool operator != (const Square_root_2_field& lh, const Square_root_2_field& rh) +{ + return !(lh == rh); +} + +template +Square_root_2_field operator * (const Int& val, const Square_root_2_field& rh) +{ + Square_root_2_field temp = rh; + temp.l *= val; + temp.r *= val; + + return temp; +} + +template +ostream& operator<<(ostream& os, const Square_root_2_field& nb) +{ + os << nb.l << " " << nb.r; + return os; +} + +// begin MT - added to please libc++ +template +bool isnan(const Square_root_2_field& x) +{ + return false; // MT TBD +} + +template +bool isinf(const Square_root_2_field& x) +{ + return false; // MT TBD +} + +template +Square_root_2_field copysign(const Square_root_2_field& x, const Square_root_2_field& y) +{ + if (y.abs()==y) + { return x;} + return -x; +} +// end MT - added to please libc++ + +#endif // CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQUARE_ROOT_2_FIELD_H + diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Traits_with_hyperbolic_offsets_adaptor.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Traits_with_hyperbolic_offsets_adaptor.h new file mode 100644 index 00000000000..92e560c13fd --- /dev/null +++ b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Traits_with_hyperbolic_offsets_adaptor.h @@ -0,0 +1,96 @@ +// Copyright (c) 1999-2004,2006-2009,2014-2015 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Monique Teillaud +// Sylvain Pion +// Andreas Fabri +// Nico Kruithof +// Manuel Caroli +// Aymeric Pellé +#ifndef CGAL_TRAITS_WITH_HYPERBOLIC_OFFSETS_ADAPTOR_H +#define CGAL_TRAITS_WITH_HYPERBOLIC_OFFSETS_ADAPTOR_H + +#include +#include + + +namespace CGAL { + +template < class K, class Functor_ > + class Traits_with_hyperbolic_offsets_adaptor { + typedef K Kernel; + typedef Functor_ Functor; + + typedef typename Kernel::Point_3 Point; + typedef HyperbolicOctagonGroup Offset; + +public: + typedef typename Kernel::Circle_2 Circle_2; + typedef typename Kernel::Construct_hyperbolic_point_2 Construct_hyperbolic_point_2; + typedef typename Functor::result_type result_type; + + Traits_with_hyperbolic_offsets_adaptor(const Circle_2 * dom) : _domain(dom) { } + + result_type operator()(const Point& p0, const Point& p1, + const Offset& o0, const Offset& o1) const { + return Functor()(pp(p0,o0),pp(p1,o1)); + } + result_type operator()(const Point& p0, const Point& p1, const Point& p2, + const Offset& o0, const Offset& o1, const Offset& o2) const { + return Functor()(pp(p0,o0),pp(p1,o1),pp(p2,o2)); + } + result_type operator()(const Point& p0, const Point& p1, + const Point& p2, const Point& p3, + const Offset& o0, const Offset& o1, + const Offset& o2, const Offset& o3) const { + return Functor()(pp(p0,o0),pp(p1,o1),pp(p2,o2),pp(p3,o3)); + } + result_type operator()(const Point& p0, const Point& p1, + const Point& p2, const Point& p3, const Point& p4, + const Offset& o0, const Offset& o1, const Offset& o2, + const Offset& o3, const Offset& o4) const { + return Functor()(pp(p0,o0),pp(p1,o1),pp(p2,o2), + pp(p3,o3),pp(p4,o4)); + } + + result_type operator()(const Point& p0, const Point& p1) const { + return Functor()(p0, p1); + } + result_type operator()(const Point& p0, const Point& p1, + const Point& p2) const { + return Functor()(p0, p1, p2); + } + result_type operator()(const Point& p0, const Point& p1, + const Point& p2, const Point& p3) const { + return Functor()(p0, p1, p2, p3); + } + result_type operator()(const Point& p0, const Point& p1, + const Point& p2, const Point& p3, const Point& p4) const { + return Functor()(p0, p1, p2, p3, p4); + } + +protected: + Point pp(const Point &p, const Offset &o) const { + return Construct_point_3(*_domain)(p,o); + } + public: + const Circle_2* _domain; +}; +} // namespace CGAL + +#endif /* CGAL_TRAITS_WITH_HYPERBOLIC_OFFSETS_ADAPTOR_H */ diff --git a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h b/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h deleted file mode 100644 index f7018f5f470..00000000000 --- a/Periodic_4_hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h +++ /dev/null @@ -1,519 +0,0 @@ -// Copyright (c) 2010 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org); you may redistribute it under -// the terms of the Q Public License version 1.0. -// See the file LICENSE.QPL distributed with CGAL. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL: svn+ssh://mbogdanov@scm.gforge.inria.fr/svn/cgal/trunk/Triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h $ -// $Id: Triangulation_hyperbolic_traits_2.h 57323 2010-07-05 10:07:39Z sloriot $ -// -// -// Author(s) : Mikhail Bogdanov - -#ifndef CGAL_TRIANGULATION_HYPERBOLIC_TRAITS_2_H -#define CGAL_TRIANGULATION_HYPERBOLIC_TRAITS_2_H - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include "boost/tuple/tuple.hpp" -#include "boost/variant.hpp" - -namespace CGAL { - -template < class R > -class Triangulation_hyperbolic_traits_2 { -public: - typedef Triangulation_hyperbolic_traits_2 Self; - - typedef R Kernel; - - typedef R Rep; - typedef typename R::RT RT; - typedef typename R::Point_2 Point_2; - typedef typename R::Vector_2 Vector_2; - typedef typename R::Triangle_2 Triangle_2; - typedef typename R::Line_2 Line_2; - typedef typename R::Ray_2 Ray_2; - - typedef typename R::Vector_3 Vector_3; - typedef typename R::Point_3 Point_3; - - typedef typename R::Less_x_2 Less_x_2; - typedef typename R::Less_y_2 Less_y_2; - typedef typename R::Compare_x_2 Compare_x_2; - typedef typename R::Compare_y_2 Compare_y_2; - typedef typename R::Orientation_2 Orientation_2; - typedef typename R::Side_of_oriented_circle_2 Side_of_oriented_circle_2; - typedef typename R::Construct_bisector_2 Construct_bisector_2; - typedef typename R::Compare_distance_2 Compare_distance_2; - typedef typename R::Construct_triangle_2 Construct_triangle_2; - typedef typename R::Construct_direction_2 Construct_direction_2; - - typedef typename R::Angle_2 Angle_2; - typedef typename R::Construct_midpoint_2 Construct_midpoint_2; - typedef typename R::Compute_squared_distance_2 Compute_squared_distance_2; - - typedef typename R::Iso_rectangle_2 Iso_rectangle_2; - typedef typename R::Circle_2 Circle_2; - - typedef boost::tuple Arc_2; - typedef typename R::Segment_2 Line_segment_2; - typedef boost::variant Segment_2; - - typedef typename R::Line_2 Euclidean_line_2; - -private: - // Poincaré disk - const Circle_2 _unit_circle; - -public: - const Circle_2& unit_circle() const - { - return _unit_circle; - } - - Angle_2 - angle_2_object() const - { return Angle_2(); } - - Compute_squared_distance_2 - compute_squared_distance_2_object() const - { return Compute_squared_distance_2(); } - - class Construct_segment_2 - { - typedef typename CGAL::Regular_triangulation_filtered_traits_2 Regular_geometric_traits_2; - typedef typename Regular_geometric_traits_2::Construct_weighted_circumcenter_2 Construct_weighted_circumcenter_2; - typedef typename Regular_geometric_traits_2::Weighted_point_2 Weighted_point_2; - typedef typename Regular_geometric_traits_2::Bare_point Bare_point; - - public: - Construct_segment_2(const Circle_2& c) : _unit_circle(c) - { - } - - Segment_2 operator()(const Point_2& p, const Point_2& q) const - { - typedef typename R::Collinear_2 Collinear_2; - if(Collinear_2()(p, q, _unit_circle.center())){ - return Line_segment_2(p, q); - } - - Weighted_point_2 wp(p); - Weighted_point_2 wq(q); - Weighted_point_2 wo(_unit_circle.center(), _unit_circle.squared_radius()); - - Bare_point center = Construct_weighted_circumcenter_2()(wp, wo, wq); - FT radius = Compute_squared_distance_2()(p, center); - - Circle_2 circle( center, radius); - // uncomment!!! - //assert(circle.has_on_boundary(p) && circle.has_on_boundary(q)); - - if(Orientation_2()(p, q, center) == LEFT_TURN) { - return Arc_2(circle, p, q); - } - return Arc_2(circle, q, p); - } - - private: - const Circle_2& _unit_circle; - }; - - Construct_segment_2 - construct_segment_2_object() const - { - return Construct_segment_2(_unit_circle); - } - - class Construct_circumcenter_2 - { - public: - Construct_circumcenter_2(const Circle_2& c) : _unit_circle(c) - {} - - // TODO: improve this function - Point_2 operator()(Point_2 p, Point_2 q, Point_2 r) - { - assert(_unit_circle.bounded_side(p) == ON_BOUNDED_SIDE); - assert(_unit_circle.bounded_side(q) == ON_BOUNDED_SIDE); - assert(_unit_circle.bounded_side(r) == ON_BOUNDED_SIDE); - - Circle_2 circle(p, q, r); - // circle must be inside the unit one - assert(CGAL::do_intersect(_unit_circle, circle) == false); - - if(circle.center() <= _unit_circle.center() && circle.center() >= _unit_circle.center()){ - return _unit_circle.center(); - } - - FT x0 = circle.center().x(), y0 = circle.center().y(); - // a*alphaˆ2 + b*alpha + c = 0; - FT a = x0*x0 + y0*y0; - FT b = a - circle.squared_radius() + _unit_circle.squared_radius(); - FT c = _unit_circle.squared_radius(); - FT D = b*b - 4*a*c; - - FT alpha = (b - CGAL::sqrt(to_double(D)))/(2*a); - - Point_2 center(x0*alpha, y0*alpha); - if(!circle.has_on_bounded_side(center)) - { std::cout << "Center does not belong to the pencil of spheres!!!" << std::endl;} ; - return center; - } - - private: - const Circle_2 _unit_circle; - }; - - Construct_circumcenter_2 - construct_circumcenter_2_object() - { - Construct_circumcenter_2(_unit_circle); - } - - Construct_midpoint_2 - construct_midpoint_2_object() const - { return Construct_midpoint_2(); } - - //for natural_neighbor_coordinates_2 - typedef typename R::FT FT; - typedef typename R::Equal_x_2 Equal_x_2; - typedef typename R::Compute_area_2 Compute_area_2; - Compute_area_2 compute_area_2_object () const - { - return Compute_area_2(); - } - - // for compatibility with previous versions - typedef Point_2 Point; - typedef Segment_2 Segment; - typedef Triangle_2 Triangle; - typedef Ray_2 Ray; - //typedef Line_2 Line; - - Triangulation_hyperbolic_traits_2() : - _unit_circle(Point_2(0, 0), 1*1) - {} - - Triangulation_hyperbolic_traits_2(FT r) : - _unit_circle(Point_2(0, 0), r*r) - {} - - Triangulation_hyperbolic_traits_2(const Triangulation_hyperbolic_traits_2 & other) : - _unit_circle(other._unit_circle) - {} - - Triangulation_hyperbolic_traits_2 &operator= - (const Triangulation_hyperbolic_traits_2 &) - { - return *this; - } - - Less_x_2 - less_x_2_object() const - { return Less_x_2();} - - Less_y_2 - less_y_2_object() const - { return Less_y_2();} - - Compare_x_2 - compare_x_2_object() const - { return Compare_x_2();} - - Compare_y_2 - compare_y_2_object() const - { return Compare_y_2();} - - Orientation_2 - orientation_2_object() const - { return Orientation_2();} - - Side_of_oriented_circle_2 - side_of_oriented_circle_2_object() const - {return Side_of_oriented_circle_2();} - - Construct_circumcenter_2 - construct_circumcenter_2_object() const - { - return Construct_circumcenter_2(_unit_circle); - } - - class Construct_hyperbolic_bisector_2 - { - public: - Construct_hyperbolic_bisector_2(const Circle_2& unit_circle) : - _unit_circle(unit_circle) {} - - Segment_2 operator()(Point_2 p, Point_2 q) const - { - // If two points are almost of the same distance to the origin, then - // the bisector is supported by the circle of huge radius etc. - // This circle is computed inexactly. - // At present time, in this case the bisector is supported by the line. - - Compute_squared_distance_2 dist = Compute_squared_distance_2(); - Point origin = _unit_circle.center(); - FT dif = dist(origin, p) - dist(origin, q); - FT eps = 0.0000000001; - - // Bisector is straight in euclidean sense - if(dif > -eps && dif < eps){ - - // ideally - //if(Compare_distance_2()(_unit_circle.center(), p, q) == EQUAL){ - - // TODO: calling R::Construct_bisector - Euclidean_line_2 l = Construct_bisector_2()(p, q); - // compute the ending points - std::pair points = find_intersection(l); - // TODO: improve - Vector_2 v(points.first, points.second); - if(v*l.to_vector() > 0){ - return Line_segment_2(points.first, points.second); - } - return Line_segment_2(points.second, points.first); - } - - Circle_2 c = construct_supporting_circle(p, q); - // compute the ending points - std::pair points = find_intersection(c); - - if(Orientation_2()(points.first, points.second, c.center()) == LEFT_TURN) { - return Arc_2(c, points.first, points.second); - } - return Arc_2(c, points.second, points.first); - } - - private: - // The cirle belongs to the pencil with limit points p and q - Circle_2 construct_supporting_circle(Point_2 p, Point_2 q) const - { - // p, q are zero-circles - // (x, y, xˆ2 + yˆ2 - rˆ2) = alpha*(xp, yp, xpˆ2 + ypˆ2) + (1-alpha)*(xq, yq, xqˆ2 + yqˆ2) - // xˆ2 + yˆ2 - rˆ2 = Rˆ2, where R - is a radius of the given unit circle - FT op = p.x()*p.x() + p.y()*p.y(); - FT oq = q.x()*q.x() + q.y()*q.y(); - FT alpha = (_unit_circle.squared_radius() - oq) / (op - oq); - - FT x = alpha*p.x() + (1-alpha)*q.x(); - FT y = alpha*p.y() + (1-alpha)*q.y(); - FT radius = x*x + y*y - _unit_circle.squared_radius(); - - //improve - typename R::Line_2 l = typename R::Construct_bisector_2()(p, q); - Point_2 middle = Construct_midpoint_2()(p, q); - Point_2 temp = middle + l.to_vector(); - if(Orientation_2()(middle, temp, Point_2(x, y)) == ON_POSITIVE_SIDE){ - return Circle_2(Point_2(x, y), radius, CLOCKWISE); - } - - return Circle_2(Point_2(x, y), radius, COUNTERCLOCKWISE); - } - - // Find intersection of an input circle orthogonal to the Poincaré disk - // and the circle representing this disk - - // TODO: sqrt(to_double()?) - std::pair find_intersection(Circle_2& circle) const - { - FT x = circle.center().x(), y = circle.center().y(); - - // axˆ2 + 2bˆx + c = 0; - FT a = x*x + y*y; - FT b = -_unit_circle.squared_radius() * x; - FT c = _unit_circle.squared_radius()*_unit_circle.squared_radius() - _unit_circle.squared_radius()*y*y; - assert(b*b - a*c > 0); - FT D = CGAL::sqrt(to_double(b*b - a*c)); - - FT x1 = (-b - D)/a; - FT x2 = (-b + D)/a; - FT y1 = (_unit_circle.squared_radius() - x1*x)/y; - FT y2 = (_unit_circle.squared_radius() - x2*x)/y; - - return std::make_pair(Point_2(x1, y1), Point_2(x2, y2)); - } - - // Find intersection of an input line orthogonal to the Poincaré disk - // and the circle representing this disk - - // TODO: sqrt(to_double()?) - std::pair find_intersection(Euclidean_line_2& l) const - { - typedef typename R::Vector_2 Vector_2; - Vector_2 v = l.to_vector(); - - // normalize the vector - FT squared_coeff = _unit_circle.squared_radius()/v.squared_length(); - FT coeff = CGAL::sqrt(to_double(squared_coeff)); - - Point_2 p1(coeff*v.x(), coeff*v.y()); - Point_2 p2(-p1.x(), -p1.y()); - return std::make_pair(p1, p2); - } - - private: - const Circle_2 _unit_circle; - }; - - Construct_hyperbolic_bisector_2 - construct_hyperbolic_bisector_2_object() const - { return Construct_hyperbolic_bisector_2(_unit_circle);} - - Construct_bisector_2 - construct_bisector_2_object() const - {return Construct_bisector_2();} - - Compare_distance_2 - compare_distance_2_object() const - {return Compare_distance_2();} - - Construct_triangle_2 construct_triangle_2_object() const - {return Construct_triangle_2();} - - Construct_direction_2 construct_direction_2_object() const - {return Construct_direction_2();} - - class Construct_ray_2 - { - public: - Construct_ray_2(Circle_2 c) : - _unit_circle(c) {} - - Segment_2 operator()(Point_2 p, Segment_2 l) const - { - if(typename R::Segment_2* s = boost::get(&l)){ - return operator()(p, *s); - } - if(Arc_2* arc = boost::get(&l)){ - if(get<0>(*arc).orientation() == CLOCKWISE){ - get<1>(*arc) = p; - return *arc; - } - get<2>(*arc) = p; - return *arc; - } - assert(false); - return Segment_2(); - } - - Segment_2 operator()(Point_2 p, typename R::Segment_2 s) const - { - return typename R::Segment_2(p, s.target()); - } - - private: - - const Circle_2 _unit_circle; - }; - - Construct_ray_2 construct_ray_2_object() const - {return Construct_ray_2(_unit_circle);} - - // For details see the rapport RR-8146 - class Is_hyperbolic - { - public: - - bool operator() (const Point& p0, const Point& p1, const Point& p2) const - { - Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), - p1.x()*p1.x() + p1.y()*p1.y(), - p2.x()*p2.x() + p2.y()*p2.y()); - - Vector_3 v1 = Vector_3(p0.x(), p1.x(), p2.x()); - Vector_3 v2 = Vector_3(p0.y(), p1.y(), p2.y()); - Vector_3 v3 = Vector_3(FT(1), FT(1), FT(1)); - - FT dt0 = CGAL::determinant(v0, v1, v3); - FT dt1 = CGAL::determinant(v0, v2, v3); - FT dt2 = CGAL::determinant(v0 - v3, v1, v2); - - return dt0*dt0 + dt1*dt1 - dt2*dt2 < 0; - } - - bool operator() (const Point& p0, const Point& p1, const Point& p2, int& ind) const - { - if (this->operator()(p0, p1, p2) == false) { - ind = find_non_hyperbolic_edge(p0, p1, p2); - return false; - } - return true; - } - - private: - - // assume the face (p0, p1, p2) is non-hyperbolic - int find_non_hyperbolic_edge(const Point& p0, const Point& p1, const Point& p2) const - { - typedef typename R::Direction_2 Direction_2; - - Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(), - p1.x()*p1.x() + p1.y()*p1.y(), - p2.x()*p2.x() + p2.y()*p2.y()); - - Vector_3 v1 = Vector_3(p0.x(), p1.x(), p2.x()); - Vector_3 v2 = Vector_3(p0.y(), p1.y(), p2.y()); - Vector_3 v3 = Vector_3(FT(1), FT(1), FT(1)); - - FT dt0 = CGAL::determinant(v0, 2*v2, -v3); - FT dt1 = CGAL::determinant(2*v1, v0, -v3); - FT dt2 = CGAL::determinant(2*v1, 2*v2, -v3); - - Direction_2 d0(p0.x()*dt2 - dt0, p0.y()*dt2 - dt1); - Direction_2 d1(p1.x()*dt2 - dt0, p1.y()*dt2 - dt1); - Direction_2 d2(p2.x()*dt2 - dt0, p2.y()*dt2 - dt1); - - Direction_2 d(dt0, dt1); - - if(d.counterclockwise_in_between(d0, d1)) { - return 2; - } - - if(d.counterclockwise_in_between(d1, d2)) { - return 0; - } - - return 1; - } - }; - - Is_hyperbolic Is_hyperbolic_object() const - { return Is_hyperbolic(); } -}; - -// Take out the code below to some separate file - -#ifdef CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H -template <> -struct Triangulation_structural_filtering_traits< Triangulation_hyperbolic_traits_2 > { - typedef Tag_true Use_structural_filtering_tag; -}; -#endif // CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H - -#ifdef CGAL_EXACT_PREDICATES_INEXACT_CONSTRUCTIONS_KERNEL_H -template <> -struct Triangulation_structural_filtering_traits< Triangulation_hyperbolic_traits_2 > { - typedef Tag_true Use_structural_filtering_tag; -}; -#endif // CGAL_EXACT_PREDICATES_INEXACT_CONSTRUCTIONS_KERNEL_H - -} //namespace CGAL - -#endif // CGAL_TRIANGULATION_HYPERBOLIC_TRAITS_2_H