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

Reorganisation of include/ directory and modifications regarding class hierarchy. Ported a periodic triangulation traits class to the hyperbolic case.

This commit is contained in:
Iordan Iordanov 2016-03-09 16:56:29 +01:00
parent e8f984caf0
commit 5996dd660b
23 changed files with 2018 additions and 2533 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/CMakeLists.txt
View File

@ -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} )

779
Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_Delaunay_triangulation_2_demo.cpp Normal file
View File

@ -0,0 +1,779 @@
#include <fstream>
// CGAL headers
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Delaunay_hyperbolic_triangulation_2.h>
#include <CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h>
// to be deleted
#include <CGAL/Qt/HyperbolicPainterOstream.h>
//
#include <CGAL/point_generators_2.h>
// Qt headers
#include <QtGui>
#include <QString>
#include <QActionGroup>
#include <QFileDialog>
#include <QInputDialog>
#include <QGraphicsEllipseItem>
#include <string>
// 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 <CGAL/Qt/TriangulationGraphicsItem.h>
#include <CGAL/Qt/VoronoiGraphicsItem.h>
// store color
#include <CGAL/TranslationInfo.h>
// visualise color
#include <CGAL/Qt/TriangulationGraphicsItemWithColorInfo.h>
// unique words
//#include <temp.h>
#include <CGAL/Square_root_2_field.h>
#include <CGAL/Hyperbolic_octagon_group.h>
#include <CGAL/Hyperbolic_random_points_in_disc_2.h>
// dummy points
//#include <CGAL/Periodic_2_hyperbolic_triangulation_dummy.h>
// for filtering
#include <set>
// for viewportsBbox
#include <CGAL/Qt/utility.h>
// the two base classes
#include "ui_Delaunay_triangulation_2.h"
#include <CGAL/Qt/DemosMainWindow.h>
#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<R> K;
typedef K::Point_2 Point;
// keep color
typedef TranslationInfo<std::string> Vb_info;
typedef CGAL::Triangulation_vertex_base_with_info_2< Vb_info, K >
Vb;
typedef CGAL::Triangulation_face_base_with_info_2 <CGAL::Hyperbolic_face_info_2, K >
Fb;
typedef CGAL::Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2<Vb, Fb> >
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<HyperbolicOctagonGroup>& v, vector<HyperbolicOctagonGroup> g, int depth = 1) {
if (depth > 1) {
recurr(v, g, depth-1);
vector<HyperbolicOctagonGroup> tmp;
vector<string> 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<Point>& p, Point input, int length) {
std::vector<HyperbolicOctagonGroup> g;
get_generators(g);
std::vector<HyperbolicOctagonGroup> v;
recurr(v, g, length);
std::set<HyperbolicOctagonGroup> 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<HyperbolicOctagonGroup>::iterator it = s.begin(); it != s.end(); it++) {
HyperbolicOctagonGroup m = *it;
pair<double, double> 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<Point>& 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<HyperbolicOctagonGroup> 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<double, double> res;
for(size_t i = 0; i < unique_words.size(); i++) {
pair<double, double> 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<Point>& points, Point input = Point(0, 0), double threshold = 6/*13.5*/, int word_length = 6/*20*/)
{
static vector<HyperbolicOctagonGroup> 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 << "," <<endl;
}
fwords << unique_words[unique_words.size() - 1].label << endl;
fwords << "];\nreturn words;";
fwords.close();
*/
//
}
ifstream findices("indicesG8_1");
long index;
vector<long> indices;
while(findices >> index) {
indices.push_back(index);
}
cout << "indices size " << indices.size() << endl;
pair<double, double> res;
for(size_t i = 0; i < indices.size(); i++) {
pair<double, double> 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<int> ccol;
Delaunay dt;
QGraphicsEllipseItem * disk;
QGraphicsScene scene;
CGAL::Qt::TriangulationGraphicsItem<Delaunay> * dgi;
CGAL::Qt::VoronoiGraphicsItem<Delaunay> * vgi;
// for drawing Voronoi diagram of the orbit of the origin
CGAL::Qt::VoronoiGraphicsItem<Delaunay> * origin_vgi;
CGAL::Qt::TriangulationMovingPoint<Delaunay> * mp;
CGAL::Qt::TriangulationConflictZone<Delaunay> * cz;
CGAL::Qt::TriangulationRemoveVertex<Delaunay> * trv;
CGAL::Qt::TriangulationPointInputAndConflictZone<Delaunay> * pi;
CGAL::Qt::TriangulationCircumcircle<Delaunay> * 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<Point> 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<Delaunay>(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<Delaunay>(&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<Delaunay>(&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<Delaunay>(&scene, &dt, this );
QObject::connect(pi, SIGNAL(generate(CGAL::Object)),
this, SLOT(processInput(CGAL::Object)));
mp = new CGAL::Qt::TriangulationMovingPoint<Delaunay>(&dt, this);
// TriangulationMovingPoint<Delaunay> 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<Delaunay>(&dt, this);
QObject::connect(trv, SIGNAL(modelChanged()),
this, SIGNAL(changed()));
tcc = new CGAL::Qt::TriangulationCircumcircle<Delaunay>(&scene, &dt, this);
tcc->setPen(QPen(Qt::red, 0, Qt::SolidLine, Qt::RoundCap, Qt::RoundJoin));
cz = new CGAL::Qt::TriangulationConflictZone<Delaunay>(&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<Point> 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<Point> 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<Point> 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_<action_name>_<signal_name>"
*/
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<K> convert;
//Circle_2 isor = convert(rect);
//CGAL::Random_points_in_disc_2<Point> 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<int>::max(),
1,
&ok);
if(!ok) {
return;
}
// wait cursor
QApplication::setOverrideCursor(Qt::WaitCursor);
//std::vector<Point> 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<Point_2> pts;
Hyperbolic_random_points_in_disc_2<GT>(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<K::Point> 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();
}

110
Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_Dirichlet_region_demo.cpp → Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_Dirichlet_region/Periodic_4_hyperbolic_triangulation_2_demo.cpp
View File

@ -6,9 +6,7 @@
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Delaunay_hyperbolic_triangulation_2.h>
/* #include <CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h> */
#include <CGAL/Triangulation_hyperbolic_traits_2.h>
#include <CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h>
// to be deleted
#include <CGAL/Qt/HyperbolicPainterOstream.h>
@ -43,8 +41,8 @@
// unique words
//#include <temp.h>
#include <CGAL/Sqrt_field.h>
#include <CGAL/Octagon_matrix.h>
#include <CGAL/Square_root_2_field.h>
#include <CGAL/Hyperbolic_octagon_group.h>
#include <CGAL/Hyperbolic_random_points_in_disc_2.h>
@ -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<R> K;
typedef K::Point_2 Point;
// keep color
typedef TranslationInfo<std::string> Vb_info;
typedef CGAL::Triangulation_vertex_base_with_info_2< Vb_info, K >
Vb;
typedef CGAL::Triangulation_face_base_with_info_2 <CGAL::Hyperbolic_face_info_2, K >
Fb;
typedef CGAL::Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2<Vb, Fb> >
/* typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2< K, CGAL::Triangulation_data_structure_2<Vb, Fb> > */
Delaunay;
typedef Delaunay::Vertex_handle Vertex_handle;
//typedef CGAL::Delaunay_hyperbolic_triangulation_2<K> Delaunay;
struct PointsComparator {
static double eps;
@ -105,12 +98,64 @@ struct PointsComparator {
double PointsComparator::eps = 0.0001;
void apply_unique_words(std::vector<Point>& 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<HyperbolicOctagonGroup>& v, vector<HyperbolicOctagonGroup> g, int depth = 1) {
if (depth > 1) {
recurr(v, g, depth-1);
vector<HyperbolicOctagonGroup> tmp;
vector<string> 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<Point>& p, Point input, int length) {
std::vector<HyperbolicOctagonGroup> g;
get_generators(g);
std::vector<HyperbolicOctagonGroup> v;
recurr(v, g, length);
std::set<HyperbolicOctagonGroup> 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<HyperbolicOctagonGroup>::iterator it = s.begin(); it != s.end(); it++) {
HyperbolicOctagonGroup m = *it;
pair<double, double> 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<Point>& 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<OctagonMatrix> unique_words;
static vector<HyperbolicOctagonGroup> 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<Point>& points, Point input = Point(0, 0), d
void apply_unique_words_G(std::vector<Point>& points, Point input = Point(0, 0), double threshold = 6/*13.5*/, int word_length = 6/*20*/)
{
static vector<OctagonMatrix> unique_words;
static vector<HyperbolicOctagonGroup> 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<Point> 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<Point> 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)
{

2
Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/CMakeLists.txt
View File

@ -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 )

9
Periodic_4_hyperbolic_triangulation_2/demo/Periodic_4_hyperbolic_triangulation_2/Periodic_4_hyperbolic_triangulation_2_demo.cpp
View File

@ -3,8 +3,8 @@
// CGAL headers
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
//#include <CGAL/Delaunay_hyperbolic_triangulation_2.h>
#include <CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h>
#include <CGAL/Delaunay_hyperbolic_triangulation_2.h>
//#include <CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h>
#include <CGAL/Triangulation_hyperbolic_traits_2.h>
// to be deleted
@ -43,7 +43,8 @@ typedef CGAL::Triangulation_hyperbolic_traits_2<R> K;
typedef K::Point_2 Point_2;
typedef K::Iso_rectangle_2 Iso_rectangle_2;
typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2<K> Delaunay;
//typedef CGAL::Periodic_4_Delaunay_hyperbolic_triangulation_2<K> Delaunay;
typedef CGAL::Delaunay_hyperbolic_triangulation_2<K> 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);
}

58
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Base_hyperbolic_translations.h
View File

@ -1,58 +0,0 @@
#ifndef CGAL_HYPERBOLIC_BASE_TRANSLATIONS_2_H
#define CGAL_HYPERBOLIC_BASE_TRANSLATIONS_2_H
#include <CGAL/Hyperbolic_isometry_2.h>
template<typename Translation>
struct Element
{
typedef typename std::list<Element> List;
typedef CGAL::Circulator_from_container<List> Circulator;
Translation g;
// circulator iterator to an inverse translation in the list
Circulator inverse;
};
template<typename Gt>
class Base_hyperbolic_translations
{
public:
typedef typename Gt::FT FT;
typedef typename std::complex<FT> complex;
typedef CGAL::Hyperbolic_isometry_2<Gt> Hyperbolic_isometry;
//typedef Element<Hyperbolic_isometry> Element_t;
//typedef std::list<Element_t> List;
//typedef typename List::iterator List_iterator;
//typedef CGAL::Circulator_from_container<List> Circulator;
typedef std::pair<Hyperbolic_isometry, int> Node;
typedef std::vector<Node> 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

125
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Cross_translations.h
View File

@ -1,125 +0,0 @@
#ifndef CGAL_HYPERBOLIC_CROSS_TRANSLATIONS_2_H
#define CGAL_HYPERBOLIC_CROSS_TRANSLATIONS_2_H
#include <CGAL/Hyperbolic_isometry_2.h>
template<typename Gt>
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<FT> complex;
typedef CGAL::Hyperbolic_isometry_2<Gt> Hyperbolic_isometry;
typedef std::vector<Hyperbolic_isometry> 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<FT> m(k1, k1);
std::complex<FT> 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

1
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Delaunay_hyperbolic_triangulation_2.h
View File

@ -134,7 +134,6 @@ public:
: Delaunay_triangulation_2<Gt,Tds>(tr)
{ CGAL_triangulation_postcondition( this->is_valid() );}
void mark_star(Vertex_handle v) const
{
if(!is_star_bounded(v)) {

132
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Diametric_translations.h
View File

@ -1,132 +0,0 @@
#ifndef CGAL_HYPERBOLIC_DIAMETRIC_TRANSLATIONS_2_H
#define CGAL_HYPERBOLIC_DIAMETRIC_TRANSLATIONS_2_H
#include <CGAL/Hyperbolic_isometry_2.h>
template<typename Gt>
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<FT> complex;
typedef CGAL::Hyperbolic_isometry_2<Gt> Hyperbolic_isometry;
typedef std::vector<Hyperbolic_isometry> 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<FT> 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<FT> > 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<FT>( 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

0
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_Delaunay_triangulation_2.h
View File

105
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_isometry_2.h
View File

@ -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 <complex>
#include <CGAL/number_utils.h>
namespace CGAL {
/*
z -> (m*z + n)/(conj(n)*z + conj(m))
*/
template<class Gt>
class Hyperbolic_isometry_2
{
public:
typedef Gt Geom_traits;
typedef typename Gt::FT FT;
typedef typename Gt::Point Point;
typedef std::complex<FT> 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 <class Gt>
std::ostream&
operator<<(std::ostream& os, const Hyperbolic_isometry_2<Gt> &isometry)
{
isometry.file_output(os);
return os;
}
} // namespace CGAL
#endif // CGAL_HYPERBOLIC_ISOMETRY_2_H

575
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Hyperbolic_octagon_group.h Normal file
View File

@ -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 <complex>
#include <iostream>
#include <vector>
#include <iterator>
#include <map>
#include <set>
#include <assert.h>
#include <string>
#include <fstream>
template<class Square_root_2_field>
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<Square_root_2_field> Self;
typedef complex<Square_root_2_field> 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<double> toComplexDouble(Matrix_element M) //const
{
Square_root_2_field rl = real(M);
Square_root_2_field img = imag(M);
return complex<double>(rl.l + sqrt(2.)*rl.r, img.l + sqrt(2.)*img.r);
}
pair<double, double> apply(double x, double y)
{
typedef complex<double> 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<double, double>(real(res), imag(res));
}
};
template<class Square_root_2_field>
Square_root_2_field Hyperbolic_octagon_group<Square_root_2_field>::factor = Square_root_2_field(-1, 1);
// just to give an order(ing)
template<class Int>
bool operator < (const complex<Square_root_2_field<Int> >& lh,
const complex<Square_root_2_field<Int> >& 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<class Square_root_2_field>
bool operator < (const Hyperbolic_octagon_group<Square_root_2_field>& lh,
const Hyperbolic_octagon_group<Square_root_2_field>& rh)
{
if (lh.A < rh.A) {
return true;
}
if (lh.A == rh.A ) {
if (lh.B < rh.B) {
return true;
}
}
return false;
}
template<class Square_root_2_field>
bool operator == (const Hyperbolic_octagon_group<Square_root_2_field>& lh,
const Hyperbolic_octagon_group<Square_root_2_field>& rh)
{
return (lh.A == rh.A && lh.B == rh.B);
}
template<class Square_root_2_field>
ostream& operator<<(ostream& os, const Hyperbolic_octagon_group<Square_root_2_field>& m)
{
os << m.A << " " << m.B;
return os;
}
typedef long long ll;
typedef Square_root_2_field<ll> SqrtField;
typedef Hyperbolic_octagon_group<SqrtField> 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<HyperbolicOctagonGroup>& 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<Element> 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<HyperbolicOctagonGroup> gens;
bool IsCanonical(const HyperbolicOctagonGroup& m);
void generate_words( set<HyperbolicOctagonGroup>& words, vector<HyperbolicOctagonGroup>& 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<HyperbolicOctagonGroup> 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<HyperbolicOctagonGroup>& 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<HyperbolicOctagonGroup m, HyperbolicOctagonGroup Aux>, m = Aux * origin * Aux^{-1}
void dfs_with_info(const pair<HyperbolicOctagonGroup, HyperbolicOctagonGroup>& new_pair,
map <HyperbolicOctagonGroup, HyperbolicOctagonGroup>& 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<HyperbolicOctagonGroup, HyperbolicOctagonGroup>(candidate, candidate_factor), visited);
}
}
}
void dfs_with_info(const HyperbolicOctagonGroup& origin,
map<HyperbolicOctagonGroup, HyperbolicOctagonGroup>& visited)
{
assert(IsCanonical(origin));
HyperbolicOctagonGroup id = HyperbolicOctagonGroup(Element(SqrtField(1, 0), SqrtField(0, 0)),
Element(SqrtField(0, 0), SqrtField(0, 0)));
pair<HyperbolicOctagonGroup, HyperbolicOctagonGroup> 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<HyperbolicOctagonGroup> visited;
set<HyperbolicOctagonGroup>::iterator it, it2;
map<long, long> nb_map;
map<long, long>::iterator mit;
dfs(m, visited);
set<pair< HyperbolicOctagonGroup, HyperbolicOctagonGroup> > 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<long, long>(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<pair< HyperbolicOctagonGroup, HyperbolicOctagonGroup> >& visited) const
{
typedef pair<HyperbolicOctagonGroup, HyperbolicOctagonGroup> 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<double> d = m.toComplexDouble(D);
complex<double> t1 = m.toComplexDouble(T1);
complex<double> t2 = m.toComplexDouble(T2);
complex<double> au = complex<double>(m.factor.l + sqrt(2.)*m.factor.r, 0);
complex<double> z1 = (t1 + sqrt(d))/(t2*sqrt(au));
complex<double> 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<HyperbolicOctagonGroup>& canonical_set, vector<HyperbolicOctagonGroup>& output)
{
set<HyperbolicOctagonGroup> redundant;
set<HyperbolicOctagonGroup>::iterator it;
for(it = canonical_set.begin(); it != canonical_set.end(); ++it) {
if(redundant.find(*it) != redundant.end()) {
continue;
}
set<HyperbolicOctagonGroup> visited;
dfs(*it, visited);
visited.erase(*it);
redundant.insert(visited.begin(), visited.end());
output.push_back(*it);
}
}
void generate_unique_words(vector<HyperbolicOctagonGroup>& output, double threshold = 10, int word_length = 13)
{
get_generators(gens);
set<HyperbolicOctagonGroup> unique_words;
vector<HyperbolicOctagonGroup> temp;
generate_words(unique_words, temp, word_length, threshold);
double l = 0;
set<HyperbolicOctagonGroup>::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<HyperbolicOctagonGroup>& 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

568
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Octagon_matrix.h
View File

@ -1,568 +0,0 @@
#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_OCTAGON_MATRIX_H
#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_OCTAGON_MATRIX_H
#include <complex>
#include <iostream>
#include <vector>
#include <iterator>
#include <map>
#include <set>
#include <assert.h>
//#include <unordered_set>
#include <string>
#include <fstream>
template<class Sqrt_field>
class Octagon_matrix
{
public:
typedef Octagon_matrix<Sqrt_field> Self;
typedef complex<Sqrt_field> 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<double> toComplexDouble(Extended_field M) //const
{
Sqrt_field rl = real(M);
Sqrt_field img = imag(M);
return complex<double>(rl.l + sqrt(2.)*rl.r, img.l + sqrt(2.)*img.r);
}
pair<double, double> apply(double x, double y)
{
typedef complex<double> 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<double, double>(real(res), imag(res));
}
//private:
static Sqrt_field factor;
};
template<class Sqrt_field>
Sqrt_field Octagon_matrix<Sqrt_field>::factor = Sqrt_field(-1, 1);
// just to give an order(ing)
template<class Int>
bool operator < (const complex<Sqrt_field<Int> >& lh,
const complex<Sqrt_field<Int> >& 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<class Sqrt_field>
bool operator < (const Octagon_matrix<Sqrt_field>& lh,
const Octagon_matrix<Sqrt_field>& rh)
{
if (lh.M11 < rh.M11) {
return true;
}
if (lh.M11 == rh.M11 ) {
if (lh.M12 < rh.M12) {
return true;
}
}
return false;
}
template<class Sqrt_field>
bool operator == (const Octagon_matrix<Sqrt_field>& lh,
const Octagon_matrix<Sqrt_field>& rh)
{
return (lh.M11 == rh.M11 && lh.M12 == rh.M12);
}
template<class Sqrt_field>
ostream& operator<<(ostream& os, const Octagon_matrix<Sqrt_field>& m)
{
os << m.M11 << " " << m.M12;
return os;
}
typedef long long ll;
typedef Sqrt_field<ll> SqrtField;
typedef Octagon_matrix<SqrtField> 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<OctagonMatrix>& 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<Entry> 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<OctagonMatrix> gens;
bool IsCanonical(const OctagonMatrix& m);
void generate_words( set<OctagonMatrix>& words, vector<OctagonMatrix>& 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<OctagonMatrix> 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<OctagonMatrix>& 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<OctagonMatrix m, OctagonMatrix Aux>, m = Aux * origin * Aux^{-1}
void dfs_with_info(const pair<OctagonMatrix, OctagonMatrix>& new_pair,
map <OctagonMatrix, OctagonMatrix>& 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<OctagonMatrix, OctagonMatrix>(candidate, candidate_factor), visited);
}
}
}
void dfs_with_info(const OctagonMatrix& origin,
map<OctagonMatrix, OctagonMatrix>& visited)
{
assert(IsCanonical(origin));
OctagonMatrix id = OctagonMatrix(Entry(SqrtField(1, 0), SqrtField(0, 0)),
Entry(SqrtField(0, 0), SqrtField(0, 0)));
pair<OctagonMatrix, OctagonMatrix> 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<OctagonMatrix> visited;
set<OctagonMatrix>::iterator it, it2;
map<long, long> nb_map;
map<long, long>::iterator mit;
dfs(m, visited);
set<pair< OctagonMatrix, OctagonMatrix> > 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<long, long>(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<pair< OctagonMatrix, OctagonMatrix> >& visited) const
{
typedef pair<OctagonMatrix, OctagonMatrix> 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<double> d = m.toComplexDouble(D);
complex<double> t1 = m.toComplexDouble(T1);
complex<double> t2 = m.toComplexDouble(T2);
complex<double> au = complex<double>(m.factor.l + sqrt(2.)*m.factor.r, 0);
complex<double> z1 = (t1 + sqrt(d))/(t2*sqrt(au));
complex<double> 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<OctagonMatrix>& canonical_set, vector<OctagonMatrix>& output)
{
set<OctagonMatrix> redundant;
set<OctagonMatrix>::iterator it;
for(it = canonical_set.begin(); it != canonical_set.end(); ++it) {
if(redundant.find(*it) != redundant.end()) {
continue;
}
set<OctagonMatrix> visited;
dfs(*it, visited);
visited.erase(*it);
redundant.insert(visited.begin(), visited.end());
output.push_back(*it);
}
}
void generate_unique_words(vector<OctagonMatrix>& output, double threshold = 10, int word_length = 13)
{
get_generators(gens);
set<OctagonMatrix> unique_words;
vector<OctagonMatrix> temp;
generate_words(unique_words, temp, word_length, threshold);
double l = 0;
set<OctagonMatrix>::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<OctagonMatrix>& 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

4
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_Delaunay_hyperbolic_triangulation_2.h
View File

@ -22,9 +22,9 @@
#include <vector>
// 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 <CGAL/Periodic_4_hyperbolic_triangulation_2.h>

66
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_construct_hyperbolic_point_2.h Normal file
View File

@ -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 <Monique.Teillaud@sophia.inria.fr>
// Sylvain Pion <Sylvain.Pion@sophia.inria.fr>
// Andreas Fabri <Andreas.Fabri@sophia.inria.fr>
// Nico Kruithof <Nico.Kruithof@sophia.inria.fr>
// Manuel Caroli <Manuel.Caroli@sophia.inria.fr>
// Aymeric Pellé <Aymeric.Pelle@sophia.inria.fr>
#ifndef CGAL_PERIODIC_4_CONSTRUCT_HYPERBOLIC_POINT_2_H
#define CGAL_PERIODIC_4_CONSTRUCT_HYPERBOLIC_POINT_2_H
#include <CGAL/Square_root_2_field.h>
#include <CGAL/Hyperbolic_octagon_group.h>
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<double, double> 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

570
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_2.h.misha
View File

@ -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 <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Hyperbolic_face_info_2.h>
#include <stack>
#include <set>
namespace CGAL {
template < class Gt,
class Tds = Triangulation_data_structure_2 <
Triangulation_vertex_base_2<Gt>,
Triangulation_face_base_with_info_2<Hyperbolic_face_info_2, Gt> > >
class Delaunay_hyperbolic_triangulation_2 : public Delaunay_triangulation_2<Gt,Tds>
{
public:
typedef Delaunay_hyperbolic_triangulation_2<Gt, Tds> Self;
typedef Delaunay_triangulation_2<Gt,Tds> Base;
typedef Triangulation_face_base_with_info_2<Hyperbolic_face_info_2, Gt> 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,Tds>(gt) {}
Delaunay_hyperbolic_triangulation_2(
const Delaunay_hyperbolic_triangulation_2<Gt,Tds> &tr)
: Delaunay_triangulation_2<Gt,Tds>(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<InputIterator>::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<InputIterator>::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<Face_handle> visited_faces;
// maintain a stack to be able to backtrack
// to the most recent faces which neighbors are not visited
std::stack<Face_handle> 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<class Mark_face_test>
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<All_faces_iterator, Infinite_hyperbolic_tester>
{
typedef Filter_iterator<All_faces_iterator,
Infinite_hyperbolic_tester> 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<All_edges_iterator,
Infinite_hyperbolic_tester> 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

96
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_Delaunay_triangulation_traits_2.h Normal file
View File

@ -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 <Nico.Kruithof@sophia.inria.fr>
// Manuel Caroli <Manuel.Caroli@sophia.inria.fr>
#ifndef CGAL_PERIODIC_4_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H
#define CGAL_PERIODIC_4_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H
//#include <CGAL/basic.h>
#include <CGAL/Traits_with_hyperbolic_offsets_adaptor.h>
#include <CGAL/Periodic_4_construct_hyperbolic_point_2.h>
#include <CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h>
#include <CGAL/Hyperbolic_octagon_group.h>
namespace CGAL {
template < class THT2 >
class Periodic_4_hyperbolic_Delaunay_triangulation_traits_base_2
: public Periodic_3_triangulation_traits_base_3<THT2>
{
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<Self,
typename K::Side_of_oriented_circle_2>
Side_of_oriented_circle_2;
typedef Traits_with_offsets_adaptor<Self, typename K::Compare_distance_2>
Compare_distance_2;
typedef Traits_with_offsets_adaptor<Self,
typename K::Side_of_bounded_disc_2>
Side_of_bounded_disc_2;
// Delaunay specific constructions
typedef Traits_with_hyperbolic_offsets_adaptor<Self,
typename K::Construct_circumcenter_2>
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

182
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_offset_2.h
View File

@ -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 <Nico@nghk.nl>
#ifndef CGAL_PERIODIC_4_HYPERBOLIC_OFFSET_2_H
#define CGAL_PERIODIC_4_HYPERBOLIC_OFFSET_2_H
#include <CGAL/basic.h>
#include <CGAL/triangulation_assertions.h>
#include <CGAL/Cartesian.h>
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 <class K2>
// 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 <class K>
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

123
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Periodic_4_hyperbolic_triangulation_traits_2.h Normal file
View File

@ -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 <Nico.Kruithof@sophia.inria.fr>
// Manuel Caroli <Manuel.Caroli@sophia.inria.fr>
#ifndef CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H
#define CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H
#include <CGAL/Traits_with_hyperbolic_offsets_adaptor.h>
#include <CGAL/Triangulation_hyperbolic_traits_2.h>
#include <CGAL/Square_root_2_field.h>
#include <CGAL/Hyperbolic_octagon_group.h>
#include <CGAL/Periodic_4_construct_hyperbolic_point_2.h>
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<Self, typename K::Compare_xy_2>
Compare_xy_2;
typedef Traits_with_hyperbolic_offsets_adaptor<Self, typename K::Orientation_2>
Orientation_2;
// Triangulation constructions
typedef Periodic_4_construct_hyperbolic_point_2<Self, typename K::Construct_point_2>
Construct_hyperbolic_point_2;
typedef Traits_with_hyperbolic_offsets_adaptor<Self, typename K::Construct_segment_2>
Construct_segment_2;
typedef Traits_with_hyperbolic_offsets_adaptor<Self, typename K::Construct_triangle_2>
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 <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
namespace CGAL
{
template < typename K >
class Periodic_4_hyperbolic_Delaunay_triangulation_traits_2;
}
#endif // CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_TRAITS_2_H

217
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Sqrt_field.h
View File

@ -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 <complex>
#include <iostream>
#include <vector>
#include <iterator>
#include <map>
#include <set>
#include <assert.h>
//#include <unordered_set>
#include <string>
#include <fstream>
//#include <boost/optional.hpp>
using namespace std;
template<class Int>
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<class Int>
bool operator < (const Sqrt_field<Int>& lh, const Sqrt_field<Int>& 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<class Int>
bool operator > (const Sqrt_field<Int>& lh, const Sqrt_field<Int>& 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<class Int>
bool operator >= (const Sqrt_field<Int>& lh, const Sqrt_field<Int>& 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<class Int>
bool operator == (const Sqrt_field<Int>& lh, const Sqrt_field<Int>& rh)
{
return (lh.l == rh.l && lh.r == rh.r);
}
template<class Int>
bool operator != (const Sqrt_field<Int>& lh, const Sqrt_field<Int>& rh)
{
return !(lh == rh);
}
template<class Int>
Sqrt_field<Int> operator * (const Int& val, const Sqrt_field<Int>& rh)
{
Sqrt_field<Int> temp = rh;
temp.l *= val;
temp.r *= val;
return temp;
}
template<class Int>
ostream& operator<<(ostream& os, const Sqrt_field<Int>& nb)
{
os << nb.l << " " << nb.r;
return os;
}
// begin MT - added to please libc++
template<class Int>
bool isnan(const Sqrt_field<Int>& x)
{
return false; // MT TBD
}
template<class Int>
bool isinf(const Sqrt_field<Int>& x)
{
return false; // MT TBD
}
template<class Int>
Sqrt_field<Int> copysign(const Sqrt_field<Int>& x, const Sqrt_field<Int>& y)
{
if (y.abs()==y)
{ return x;}
return -x;
}
// end MT - added to please libc++
#endif // CGAL_PERIODIC_4_HYPERBOLIC_TRIANGULATION_SQRT_FIELD_H

199
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Square_root_2_field.h Normal file
View File

@ -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 <complex>
#include <iostream>
#include <vector>
#include <iterator>
#include <map>
#include <set>
#include <assert.h>
#include <string>
#include <fstream>
using namespace std;
template<class Int>
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<class Int>
bool operator < (const Square_root_2_field<Int>& lh, const Square_root_2_field<Int>& 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<class Int>
bool operator >= (const Square_root_2_field<Int>& lh, const Square_root_2_field<Int>& 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<class Int>
bool operator == (const Square_root_2_field<Int>& lh, const Square_root_2_field<Int>& rh)
{
return (lh.l == rh.l && lh.r == rh.r);
}
template<class Int>
bool operator != (const Square_root_2_field<Int>& lh, const Square_root_2_field<Int>& rh)
{
return !(lh == rh);
}
template<class Int>
Square_root_2_field<Int> operator * (const Int& val, const Square_root_2_field<Int>& rh)
{
Square_root_2_field<Int> temp = rh;
temp.l *= val;
temp.r *= val;
return temp;
}
template<class Int>
ostream& operator<<(ostream& os, const Square_root_2_field<Int>& nb)
{
os << nb.l << " " << nb.r;
return os;
}
// begin MT - added to please libc++
template<class Int>
bool isnan(const Square_root_2_field<Int>& x)
{
return false; // MT TBD
}
template<class Int>
bool isinf(const Square_root_2_field<Int>& x)
{
return false; // MT TBD
}
template<class Int>
Square_root_2_field<Int> copysign(const Square_root_2_field<Int>& x, const Square_root_2_field<Int>& 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

96
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Traits_with_hyperbolic_offsets_adaptor.h Normal file
View File

@ -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 <Monique.Teillaud@sophia.inria.fr>
// Sylvain Pion <Sylvain.Pion@sophia.inria.fr>
// Andreas Fabri <Andreas.Fabri@sophia.inria.fr>
// Nico Kruithof <Nico.Kruithof@sophia.inria.fr>
// Manuel Caroli <Manuel.Caroli@sophia.inria.fr>
// Aymeric Pellé <Aymeric.Pelle@sophia.inria.fr>
#ifndef CGAL_TRAITS_WITH_HYPERBOLIC_OFFSETS_ADAPTOR_H
#define CGAL_TRAITS_WITH_HYPERBOLIC_OFFSETS_ADAPTOR_H
#include <CGAL/Square_root_2_field.h>
#include <CGAL/Hyperbolic_octagon_group.h>
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 */

519
Periodic_4_hyperbolic_triangulation_2/include/CGAL/Triangulation_hyperbolic_traits_2.h
View File

@ -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 <CGAL/Point_2.h>
#include <CGAL/Segment_2.h>
#include <CGAL/Triangle_2.h>
#include <CGAL/Line_2.h>
#include <CGAL/Ray_2.h>
#include <CGAL/predicates_on_points_2.h>
#include <CGAL/basic_constructions_2.h>
#include <CGAL/distance_predicates_2.h>
#include <CGAL/Regular_triangulation_filtered_traits_2.h>
#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<R> 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<Circle_2, Point_2, Point_2> Arc_2;
typedef typename R::Segment_2 Line_segment_2;
typedef boost::variant<Arc_2, Line_segment_2> 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<R> 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<Point_2, Point_2> 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<Point_2, Point_2> 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<Point_2, Point_2> 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<Point_2, Point_2> 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<typename R::Segment_2>(&l)){
return operator()(p, *s);
}
if(Arc_2* arc = boost::get<Arc_2>(&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<Epeck> > {
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<Epick> > {
typedef Tag_true Use_structural_filtering_tag;
};
#endif // CGAL_EXACT_PREDICATES_INEXACT_CONSTRUCTIONS_KERNEL_H
} //namespace CGAL
#endif // CGAL_TRIANGULATION_HYPERBOLIC_TRAITS_2_H