From 7b2744b4c5c9ca9e5e1acb8fddefce87c1510088 Mon Sep 17 00:00:00 2001 From: POUGET Marc Date: Fri, 10 Jan 2025 21:59:44 +0100 Subject: [PATCH] test of new design with a base class traingulation and a Delaunay one --- .../Triangulation_on_hyperbolic_surface_2.cpp | 44 ++ ...ay_Triangulation_on_hyperbolic_surface_2.h | 229 +++++++ .../Triangulation_on_hyperbolic_surface_2.h | 626 ++++++++++++++++++ 3 files changed, 899 insertions(+) create mode 100644 Hyperbolic_surface_triangulation_2/examples/Hyperbolic_surface_triangulation_2/Triangulation_on_hyperbolic_surface_2.cpp create mode 100644 Hyperbolic_surface_triangulation_2/include/CGAL/Delaunay_Triangulation_on_hyperbolic_surface_2.h create mode 100644 Hyperbolic_surface_triangulation_2/include/CGAL/Triangulation_on_hyperbolic_surface_2.h diff --git a/Hyperbolic_surface_triangulation_2/examples/Hyperbolic_surface_triangulation_2/Triangulation_on_hyperbolic_surface_2.cpp b/Hyperbolic_surface_triangulation_2/examples/Hyperbolic_surface_triangulation_2/Triangulation_on_hyperbolic_surface_2.cpp new file mode 100644 index 000000000000..43ccaa96107c --- /dev/null +++ b/Hyperbolic_surface_triangulation_2/examples/Hyperbolic_surface_triangulation_2/Triangulation_on_hyperbolic_surface_2.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Exact_rational Rational; +typedef CGAL::Simple_cartesian Kernel; +typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2 ParentTraits; +typedef CGAL::Hyperbolic_surface_traits_2 Traits; +typedef CGAL::Hyperbolic_fundamental_domain_2 Domain; +typedef CGAL::Hyperbolic_fundamental_domain_factory_2 Factory; +typedef CGAL::Triangulation_on_hyperbolic_surface_2 Triangulation; +typedef CGAL::Delaunay_triangulation_on_hyperbolic_surface_2 Delaunay_triangulation; + +int main(){ + // Generates the domain: + Factory factory = Factory(); + Domain domain = factory.make_hyperbolic_fundamental_domain_g2(time(NULL)); + + // Triangulates the domain: + Triangulation triangulation = Triangulation(domain); + + // Saves the triangulation: + std::ofstream output_file = std::ofstream ("OutputTriangulation.txt"); + output_file << triangulation; + output_file.close(); + + // Prints the triangulation: + std::cout << triangulation << std::endl; + + // Generates a Delaunay triangulation + Delaunay_triangulation dt = Delaunay_triangulation(domain); + + triangulation.has_anchor(); + dt.has_anchor(); + dt.make_Delaunay(); + // Prints the triangulation: + std::cout << dt << std::endl; + + return 0; +} diff --git a/Hyperbolic_surface_triangulation_2/include/CGAL/Delaunay_Triangulation_on_hyperbolic_surface_2.h b/Hyperbolic_surface_triangulation_2/include/CGAL/Delaunay_Triangulation_on_hyperbolic_surface_2.h new file mode 100644 index 000000000000..a320293ed749 --- /dev/null +++ b/Hyperbolic_surface_triangulation_2/include/CGAL/Delaunay_Triangulation_on_hyperbolic_surface_2.h @@ -0,0 +1,229 @@ +// Copyright (c) 2024 +// INRIA Nancy (France), and Université Gustave Eiffel Marne-la-Vallee (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Vincent Despré, Loïc Dubois, Marc Pouget, Monique Teillaud + +// This file contains the declaration and the implementation of the class Delaunay_triangulation_on_hyperbolic_surface_2 + +#ifndef CGAL_DELAUNAY_HYPERBOLIC_SURFACE_TRIANGULATION_2 +#define CGAL_DELAUNAY_HYPERBOLIC_SURFACE_TRIANGULATION_2 + +#include +#include +#include + +#include +#include +#include + +namespace CGAL { + +/* +Represents a geodesic Delaunay triangulation of a closed orientable hyperbolic surface. +*/ + +template> + class Delaunay_triangulation_on_hyperbolic_surface_2: public Triangulation_on_hyperbolic_surface_2{ + public: + typedef Triangulation_on_hyperbolic_surface_2 Base;//or T_on_HS_2 + typedef typename Base::Combinatorial_map_with_cross_ratios Combinatorial_map_with_cross_ratios; + typedef typename Base::Anchor Anchor; + + typedef typename Combinatorial_map_with_cross_ratios::Dart_descriptor Dart_descriptor; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<0> Vertex_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<1> Edge_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<2> Face_range; + + typedef typename Combinatorial_map_with_cross_ratios::Dart_const_descriptor Dart_const_descriptor; + typedef typename Combinatorial_map_with_cross_ratios::Dart_const_range Dart_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<0> Vertex_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<1> Edge_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<2> Face_const_range; + + typedef typename Traits::FT Number; + typedef typename Traits::Complex Complex_number; + typedef typename Traits::Hyperbolic_point_2 Point; + typedef Hyperbolic_isometry_2 Isometry; + typedef Hyperbolic_fundamental_domain_2 Domain; + + Delaunay_triangulation_on_hyperbolic_surface_2() {}; + Delaunay_triangulation_on_hyperbolic_surface_2(const Hyperbolic_fundamental_domain_2& domain); + Delaunay_triangulation_on_hyperbolic_surface_2(Combinatorial_map_with_cross_ratios& cmap, Anchor& anchor); + + bool is_Delaunay_flippable(Dart_const_descriptor dart) const; + void flip(Dart_descriptor dart); + bool is_Delaunay() const; + int make_Delaunay(); + +protected: + Dart_descriptor pick_edge_to_flip(); + Dart_const_descriptor pick_edge_to_flip() const; +}; + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + +template +Delaunay_triangulation_on_hyperbolic_surface_2::Delaunay_triangulation_on_hyperbolic_surface_2(const Domain& domain) + : Base(domain) +{ + make_Delaunay(); +} + +template + Delaunay_triangulation_on_hyperbolic_surface_2::Delaunay_triangulation_on_hyperbolic_surface_2(Combinatorial_map_with_cross_ratios& cmap, Anchor& anchor) + : Base(cmap, anchor) + { + make_Delaunay(); + } + +//////////////////////////////////////////////////////////////////////////////// +template +bool Delaunay_triangulation_on_hyperbolic_surface_2::is_Delaunay_flippable(Dart_const_descriptor dart) const{ + CGAL_precondition(Base::is_valid()); + return ( Base::get_cross_ratio(dart).imag()>Number(0) ); +} + +template +void Delaunay_triangulation_on_hyperbolic_surface_2::flip(Dart_descriptor dart){ + CGAL_precondition(Base::is_valid()); + // First gather all the information needed + + Dart_descriptor a = Base::opposite(dart); // Get a fresh descriptor + Dart_descriptor b = Base::ccw(a); + Dart_descriptor c = Base::cw(a); + + Dart_descriptor d = Base::opposite(a); + Dart_descriptor e = Base::ccw(d); + Dart_descriptor f = Base::cw(d); + + Complex_number cross_ratio_AB = Base::get_cross_ratio(e); + Complex_number cross_ratio_BC = Base::get_cross_ratio(f); + Complex_number cross_ratio_CD = Base::get_cross_ratio(b); + Complex_number cross_ratio_DA = Base::get_cross_ratio(c); + Complex_number cross_ratio_AC = Base::get_cross_ratio(a); + + // Modify the anchor + + if (this->_anchor.dart == a){ + this->_anchor.dart = e; + this->_anchor.vertices[1] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[1], this->_anchor.vertices[2], this->_anchor.vertices[0], cross_ratio_AC)); + } else if (this->_anchor.dart == b){ + this->_anchor.vertices[2] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[0], this->_anchor.vertices[1], this->_anchor.vertices[2], cross_ratio_AC)); + } else if (this->_anchor.dart == c){ + this->_anchor.vertices[2] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[2], this->_anchor.vertices[0], this->_anchor.vertices[1], cross_ratio_AC)); + } else if (this->_anchor.dart == d){ + this->_anchor.dart = b; + this->_anchor.vertices[1] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[1], this->_anchor.vertices[2], this->_anchor.vertices[0], cross_ratio_AC)); + } else if (this->_anchor.dart == e){ + this->_anchor.vertices[2] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[0], this->_anchor.vertices[1], this->_anchor.vertices[2], cross_ratio_AC)); + } else if (this->_anchor.dart == f){ + this->_anchor.vertices[2] = Point(Base::fourth_point_from_cross_ratio(this->_anchor.vertices[2], this->_anchor.vertices[0], this->_anchor.vertices[1], cross_ratio_AC)); + } + + // Compute the new cross ratios + + Complex_number one (Number(1), Number(0)); + Complex_number cross_ratio_BD = (cross_ratio_AC) / ((cross_ratio_AC) - one) ; + Complex_number cross_ratio_AB_2 = one - (one - (cross_ratio_AB)) * (cross_ratio_AC) ; + Complex_number cross_ratio_BC_2 = one - (one - (cross_ratio_BC)) / (cross_ratio_BD) ; + Complex_number cross_ratio_CD_2 = one - (one - (cross_ratio_CD)) * (cross_ratio_AC) ; + Complex_number cross_ratio_DA_2 = one - (one - (cross_ratio_DA)) / (cross_ratio_BD) ; + + // Make the topological flip + + this->_combinatorial_map.template unlink_beta<1>(a); + this->_combinatorial_map.template unlink_beta<1>(b); + this->_combinatorial_map.template unlink_beta<1>(c); + + this->_combinatorial_map.template unlink_beta<1>(d); + this->_combinatorial_map.template unlink_beta<1>(e); + this->_combinatorial_map.template unlink_beta<1>(f); + + + this->_combinatorial_map.template link_beta<1>(b, a); + this->_combinatorial_map.template link_beta<1>(a, f); + this->_combinatorial_map.template link_beta<1>(f, b); + + this->_combinatorial_map.template link_beta<1>(e, d); + this->_combinatorial_map.template link_beta<1>(d, c); + this->_combinatorial_map.template link_beta<1>(c, e); + + // And give the new cross ratios to the edges + + this->_combinatorial_map.template info<1>(a) = cross_ratio_BD; + this->_combinatorial_map.template info<1>(e) = cross_ratio_AB_2; + this->_combinatorial_map.template info<1>(f) = cross_ratio_BC_2; + this->_combinatorial_map.template info<1>(b) = cross_ratio_CD_2; + this->_combinatorial_map.template info<1>(c) = cross_ratio_DA_2; + + // Take care of the particular cases where we need to "flip again" + + if (Base::opposite(e) == b){ + this->_combinatorial_map.template info<1>(e) = one - (one - cross_ratio_AB_2) * (cross_ratio_AC) ; + } + + if (Base::opposite(f) == c){ + this->_combinatorial_map.template info<1>(f) = one - (one - cross_ratio_BC_2) / (cross_ratio_BD) ; + } +} + +template +bool Delaunay_triangulation_on_hyperbolic_surface_2::is_Delaunay() const{ + if (! Base::is_valid()){ + return false; + } + return (pick_edge_to_flip() == nullptr); +} + +template +int Delaunay_triangulation_on_hyperbolic_surface_2::make_Delaunay(){ + CGAL_precondition(this->is_valid()); + int number_of_flips_done = 0; + + Dart_descriptor edge_to_flip = pick_edge_to_flip(); + while (edge_to_flip != nullptr){ + flip(edge_to_flip); + edge_to_flip = pick_edge_to_flip(); + number_of_flips_done++; + } + + return number_of_flips_done; +} + +//////////////////////////////////////////////////////////////////////////////// +template + typename Delaunay_triangulation_on_hyperbolic_surface_2::Dart_descriptor Delaunay_triangulation_on_hyperbolic_surface_2::pick_edge_to_flip(){ + auto &cm=this->_combinatorial_map.darts(); + for (auto it = cm.begin(); it != cm.end(); ++it){ + if ( is_Delaunay_flippable(it) ){ + return it; + } + } + return nullptr; +} + +//////////////////////////////////////////////////////////////////////////////// + +template + typename Delaunay_triangulation_on_hyperbolic_surface_2::Dart_const_descriptor Delaunay_triangulation_on_hyperbolic_surface_2::pick_edge_to_flip() const{ + const auto &cm=this->_combinatorial_map.darts(); + for (auto it = cm.begin(); it != cm.end(); ++it){ + if ( is_Delaunay_flippable(it) ){ + return it; + } + } + return nullptr; +} + +} // namespace CGAL + +#endif // CGAL_DELAUNAY_HYPERBOLIC_SURFACE_TRIANGULATION_2 diff --git a/Hyperbolic_surface_triangulation_2/include/CGAL/Triangulation_on_hyperbolic_surface_2.h b/Hyperbolic_surface_triangulation_2/include/CGAL/Triangulation_on_hyperbolic_surface_2.h new file mode 100644 index 000000000000..6f2064ba73a8 --- /dev/null +++ b/Hyperbolic_surface_triangulation_2/include/CGAL/Triangulation_on_hyperbolic_surface_2.h @@ -0,0 +1,626 @@ +// Copyright (c) 2024 +// INRIA Nancy (France), and Université Gustave Eiffel Marne-la-Vallee (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Vincent Despré, Loïc Dubois, Marc Pouget, Monique Teillaud + +// This file contains the declaration and the implementation of the class Triangulation_on_hyperbolic_surface_2 + +#ifndef CGAL_TRIANGULATION_ON_HYPERBOLIC_SURFACE_2 +#define CGAL_TRIANGULATION_ON_HYPERBOLIC_SURFACE_2 + +#include +#include +#include + +#include +#include +#include + +namespace CGAL { + +/* +Represents a geodesic triangulation of a closed orientable hyperbolic surface. +The triangulation is stored as combinatorial map decorated with one cross-ratio per edge. +It is also possible to specify an anchor for the triangulation. An anchor consists in 1) a dart of the combinatorial map, belonging by definition to a vertex V and a triangle T, together with +2) three points A,B,C in the hyperbolic plane. The points A,B,C are the three vertices in counter-clockwise order of a triangle. This triangle is a lift +of T, and A is a lift of V. +*/ + +template +struct Combinatorial_map_with_cross_ratios_item{ + template + struct Dart_wrapper{ + typedef Cell_attribute> Edge_attrib; + typedef std::tuple Attributes; + }; + }; + +template> +class Triangulation_on_hyperbolic_surface_2{ +public: + + typedef Combinatorial_map<2,Attributes> Combinatorial_map_with_cross_ratios; + + struct Anchor{ + typename Combinatorial_map_with_cross_ratios::Dart_descriptor dart; + typename Traits::Hyperbolic_point_2 vertices[3]; + }; + + typedef typename Combinatorial_map_with_cross_ratios::Dart_descriptor Dart_descriptor; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<0> Vertex_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<1> Edge_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_range<2> Face_range; + + typedef typename Combinatorial_map_with_cross_ratios::Dart_const_descriptor Dart_const_descriptor; + typedef typename Combinatorial_map_with_cross_ratios::Dart_const_range Dart_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<0> Vertex_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<1> Edge_const_range; + typedef typename Combinatorial_map_with_cross_ratios::template One_dart_per_cell_const_range<2> Face_const_range; + + typedef typename Traits::FT Number; + typedef typename Traits::Complex Complex_number; + typedef typename Traits::Hyperbolic_point_2 Point; + typedef Hyperbolic_isometry_2 Isometry; + typedef Hyperbolic_fundamental_domain_2 Domain; + + Triangulation_on_hyperbolic_surface_2() {}; + Triangulation_on_hyperbolic_surface_2(const Hyperbolic_fundamental_domain_2& domain); + Triangulation_on_hyperbolic_surface_2(Combinatorial_map_with_cross_ratios& cmap, Anchor& anchor); + + Triangulation_on_hyperbolic_surface_2& operator=(Triangulation_on_hyperbolic_surface_2 other); + + Combinatorial_map_with_cross_ratios& combinatorial_map(); + bool has_anchor() const; + Anchor& anchor(); + const Anchor& anchor() const; + + void to_stream(std::ostream& s) const; + void from_stream(std::istream& s); + + std::vector> lift(bool center=true) const; + bool is_valid() const; + + //The following methods are not documented but they are non private for internal future use. + + Dart_descriptor ccw(Dart_descriptor dart); + Dart_descriptor cw(Dart_descriptor dart); + Dart_descriptor opposite(Dart_descriptor dart); + Dart_const_descriptor const_ccw(Dart_const_descriptor dart) const; + Dart_const_descriptor const_cw(Dart_const_descriptor dart) const; + Dart_const_descriptor const_opposite(Dart_const_descriptor dart) const; + + Complex_number get_cross_ratio(Dart_const_descriptor dart) const; + + // Returns the cross ratio of the points a,b,c,d + Complex_number cross_ratio(const Point& a, const Point& b, const Point& c, const Point& d) const; + // Returns the point d such that the cross ratio of a,b,c,d is cratio + Point fourth_point_from_cross_ratio(const Point& a, const Point& b, const Point& c, const Complex_number& cratio) const; + +// Wrapper around the Cmap for iterating over vertices, edges or faces. +Vertex_range vertices_range(){ + return _combinatorial_map.template one_dart_per_cell<0>(); +} +Edge_range edges_range(){ + return _combinatorial_map.template one_dart_per_cell<1>(); +} +Face_range faces_range(){ + return _combinatorial_map.template one_dart_per_cell<2>(); +} +Vertex_const_range vertices_const_range() const { + return _combinatorial_map.template one_dart_per_cell<0>(); +} +Edge_const_range edges_const_range() const { + return _combinatorial_map.template one_dart_per_cell<1>(); +} +Face_const_range faces_const_range() const { + return _combinatorial_map.template one_dart_per_cell<2>(); +} + +protected: + Combinatorial_map_with_cross_ratios _combinatorial_map; + bool _has_anchor = false; + Anchor _anchor; + + Dart_descriptor pick_edge_to_flip(); + Dart_const_descriptor pick_edge_to_flip() const; + + void copy_from(Combinatorial_map_with_cross_ratios& cmap); + void copy_from(Combinatorial_map_with_cross_ratios& cmap, const Anchor& anchor); +}; + + template std::ostream& operator<<(std::ostream& s, const Triangulation_on_hyperbolic_surface_2& triangulation); + template void operator>>(std::istream& s, Triangulation_on_hyperbolic_surface_2& triangulation); + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// + +template +Triangulation_on_hyperbolic_surface_2::Triangulation_on_hyperbolic_surface_2(const Domain& domain){ + // (Triangulates by adding an internal edge between domain.vertex(size-1) and the other vertices) + _combinatorial_map.clear(); + int size = domain.size(); + + // Make the triangles + std::vector dart_of_triangle(size-2); + for (int k=0; k(dart_1, dart_2); + _combinatorial_map.template set_attribute<1>(dart_1, _combinatorial_map.template create_attribute<1>(cross_ratio(p0,p1,p2,p3))); + } + + // Sew the boundary edges and set their cross ratios + for (int k1=0; k1(dart_1, dart_2)){ + _combinatorial_map.template sew<2>(dart_1, dart_2); + _combinatorial_map.template set_attribute<1>(dart_1, _combinatorial_map.template create_attribute<1>(cross_ratio(p0,p1,p2,p3))); + } + } + + // Set the anchor + _anchor.dart = dart_of_triangle[0]; + _anchor.vertices[0] = domain.vertex(size-1); + _anchor.vertices[1] = domain.vertex(0); + _anchor.vertices[2] = domain.vertex(1); + _has_anchor = true; +} + +/* template */ +/* Triangulation_on_hyperbolic_surface_2::Triangulation_on_hyperbolic_surface_2(Combinatorial_map_with_cross_ratios& cmap){ */ +/* copy_from(cmap); */ +/* } */ + +template + Triangulation_on_hyperbolic_surface_2::Triangulation_on_hyperbolic_surface_2(Combinatorial_map_with_cross_ratios& cmap, Anchor& anchor){ + copy_from(cmap, anchor); +} + +//////////////////////////////////////////////////////////////////////////////// + +template + Triangulation_on_hyperbolic_surface_2& Triangulation_on_hyperbolic_surface_2::operator=(Triangulation_on_hyperbolic_surface_2 other){ + CGAL_precondition(other->is_valid()); + if (other.has_anchor()){ + copy_from(other.combinatorial_map(), other.anchor()); + } + else { + copy_from(other.combinatorial_map()); + } + return *this; +} + +//////////////////////////////////////////////////////////////////////////////// + +template +typename Triangulation_on_hyperbolic_surface_2::Combinatorial_map_with_cross_ratios& Triangulation_on_hyperbolic_surface_2::combinatorial_map(){ + return _combinatorial_map; +} + +template +bool Triangulation_on_hyperbolic_surface_2::has_anchor() const { + CGAL_precondition(is_valid()); + return _has_anchor; +} + +template +typename Triangulation_on_hyperbolic_surface_2::Anchor& +Triangulation_on_hyperbolic_surface_2::anchor() { + CGAL_precondition(is_valid() && has_anchor()); + return _anchor; +} + +template +const typename Triangulation_on_hyperbolic_surface_2::Anchor& +Triangulation_on_hyperbolic_surface_2::anchor() const { + CGAL_precondition(is_valid() && has_anchor()); + return _anchor; +} + +template +std::vector::Dart_const_descriptor, typename Triangulation_on_hyperbolic_surface_2::Point, typename Triangulation_on_hyperbolic_surface_2::Point, typename Triangulation_on_hyperbolic_surface_2::Point>> Triangulation_on_hyperbolic_surface_2::lift(bool center) const{ + CGAL_precondition(is_valid() && has_anchor()); + std::vector> realizations; + + size_t visited_darts_mark = _combinatorial_map.get_new_mark(); + _combinatorial_map.unmark_all(visited_darts_mark); + + struct Compare { + bool operator()(std::pair const & x, std::pair const & y) { + return x.second > y.second; + } + }; + std::priority_queue, std::vector>, Compare> queue; + + std::unordered_map positions; + + Dart_const_range darts = _combinatorial_map.darts(); + + _combinatorial_map.mark(_anchor.dart, visited_darts_mark); + _combinatorial_map.mark(const_ccw(_anchor.dart), visited_darts_mark); + _combinatorial_map.mark(const_cw(_anchor.dart), visited_darts_mark); + + if (center){ + Isometry center_the_drawing = hyperbolic_translation(_anchor.vertices[0]); + positions[_anchor.dart] = center_the_drawing.evaluate(_anchor.vertices[0]); + positions[const_ccw(_anchor.dart)] = center_the_drawing.evaluate(_anchor.vertices[1]); + positions[const_cw(_anchor.dart)] = center_the_drawing.evaluate(_anchor.vertices[2]); + } else { + positions[_anchor.dart] = _anchor.vertices[0]; + positions[const_ccw(_anchor.dart)] = _anchor.vertices[1]; + positions[const_cw(_anchor.dart)] = _anchor.vertices[2]; + } + + std::tuple value = std::make_tuple(_anchor.dart, positions[_anchor.dart], positions[const_ccw(_anchor.dart)], positions[const_cw(_anchor.dart)]); + realizations.push_back(value); + + Complex_number anchor_z0 (_anchor.vertices[0].x(), _anchor.vertices[0].y()); + Complex_number anchor_z1 (_anchor.vertices[1].x(), _anchor.vertices[1].y()); + Complex_number anchor_z2 (_anchor.vertices[2].x(), _anchor.vertices[2].y()); + + double weight_of_anchor_dart = CGAL::to_double(norm(anchor_z0) + norm(anchor_z1)); + double weight_of_ccw_anchor_dart = CGAL::to_double(norm(anchor_z1) + norm(anchor_z2)); + double weight_of_cw_anchor_dart = CGAL::to_double(norm(anchor_z2) + norm(anchor_z0)); + + queue.push(std::make_pair(_anchor.dart, weight_of_anchor_dart)); + queue.push(std::make_pair(const_ccw(_anchor.dart), weight_of_ccw_anchor_dart)); + queue.push(std::make_pair(const_cw(_anchor.dart), weight_of_cw_anchor_dart)); + + + + while( ! queue.empty() ){ + Dart_const_descriptor invader = queue.top().first; + queue.pop(); + + Dart_const_descriptor invaded = const_opposite(invader); + + if (!_combinatorial_map.is_marked(invaded, visited_darts_mark)){ + _combinatorial_map.mark(invaded, visited_darts_mark); + _combinatorial_map.mark(const_ccw(invaded), visited_darts_mark); + _combinatorial_map.mark(const_cw(invaded), visited_darts_mark); + + const Point &a = positions[const_ccw(invader)]; + const Point &b = positions[const_cw(invader)]; + const Point &c = positions[invader]; + Complex_number cross_ratio = get_cross_ratio(invader); + + positions[invaded] = a; + positions[const_ccw(invaded)] = c; + Point d = fourth_point_from_cross_ratio(a, b, c, cross_ratio); + positions[const_cw(invaded)] = d; + + Complex_number za (a.x(), a.y()); + Complex_number zc (c.x(), c.y()); + double invaded_distance_to_zero = CGAL::to_double(norm(za)); + double invaded_ccw_distance_to_zero = CGAL::to_double(norm(zc)); + Complex_number znew (positions[const_cw(invaded)].x(), positions[const_cw(invaded)].y()); + double invaded_cw_distance_to_zero = CGAL::to_double(norm(znew)); + + double invaded_ccw_weight = invaded_ccw_distance_to_zero + invaded_cw_distance_to_zero; + double invaded_cw_weight = invaded_cw_distance_to_zero + invaded_distance_to_zero; + + queue.push( std::make_pair(const_ccw(invaded), invaded_ccw_weight) ); + queue.push( std::make_pair(const_cw(invaded), invaded_cw_weight) ); + + value = std::make_tuple(invaded, Point(a), Point(c), Point(d)); + realizations.push_back(value); + } + } + + _combinatorial_map.free_mark(visited_darts_mark); + + return realizations; +} + +//////////////////////////////////////////////////////////////////////////////// + +template +bool Triangulation_on_hyperbolic_surface_2::is_valid() const{ + // 1. Check the combinatorial map + + // Check that the combinatorial map is valid + if ( !_combinatorial_map.is_valid() ){ + return false; + } + + // Check that the combinatorial map has no 1,2-boundary + for (int k=1; k<3; k++){ + if ( !_combinatorial_map.is_without_boundary(k) ){ + return false; + } + } + + // 2. Check the anchor, if any + + if (_has_anchor){ + // Check that the dart descriptor of the anchor points to a dart of the combinatorial map + if ( !_combinatorial_map.is_dart_used(_anchor.dart) ){ + return false; + } + + // Check that the three vertices of the anchor lie within the open unit disk + for (int k=0; k<3; k++){ + // if (_anchor.vertices[k].get_z() >= Number(1)){ + if ( norm(Complex_number(_anchor.vertices[k].x(),_anchor.vertices[k].y())) >= Number(1)){ + return false; + } + } + } + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// + +template +void Triangulation_on_hyperbolic_surface_2::to_stream(std::ostream& s) const{ + CGAL_precondition(is_valid() && has_anchor()); + // Give indices to the darts + std::map darts_indices; + int current_dart_index = 0; + for (typename Dart_const_range::const_iterator it=_combinatorial_map.darts().begin(); it!=_combinatorial_map.darts().end(); ++it){ + darts_indices[it] = current_dart_index; + current_dart_index++; + } + + // Store the number of darts + s << current_dart_index << std::endl; + + // Store the anchor, if any + if (_has_anchor){ + s << "yes" << std::endl; + s << darts_indices[_anchor.dart] << std::endl; + s << _anchor.vertices[0] << std::endl; + s << _anchor.vertices[1] << std::endl; + s << _anchor.vertices[2] << std::endl; + } else { + s << "no" << std::endl; + } + + // Store the triangles + for (typename Face_const_range::const_iterator it = faces_const_range().begin(); it != faces_const_range().end(); ++it){ + s << darts_indices[it] << std::endl; + s << darts_indices[const_cw(it)] << std::endl; + s << darts_indices[const_ccw(it)] << std::endl; + } + + // Store the edges + for (typename Edge_const_range::const_iterator it = edges_const_range().begin(); it != edges_const_range().end(); ++it){ + s << darts_indices[it] << std::endl; + s << darts_indices[const_opposite(it)] << std::endl; + s << get_cross_ratio(it); + } +} + +template +void Triangulation_on_hyperbolic_surface_2::from_stream(std::istream& s){ + _combinatorial_map.clear(); + + // Load the number of darts + std::string line; + s >> line; + int nb_darts = std::stoi(line); + + // Load the anchor + int anchor_dart_id; + s >> line; + if (!line.compare("yes")){ + _has_anchor = true; + + s >> line; + anchor_dart_id = std::stoi(line); // (*) _anchor.dart_id is set at the end of the function + + s >> _anchor.vertices[0]; + s >> _anchor.vertices[1]; + s >> _anchor.vertices[2]; + } else { + _has_anchor = false; + } + + // Load the triangles + std::vector darts_by_id (nb_darts); + int index1, index2, index3; + for (int k=0; k> line; + index1 = std::stoi(line); + s >> line; + index2 = std::stoi(line); + s >> line; + index3 = std::stoi(line); + + darts_by_id[index1] = triangle_dart; + darts_by_id[index2] = cw(triangle_dart); + darts_by_id[index3] = ccw(triangle_dart); + } + + // Load the edges + Dart_descriptor dart_1, dart_2; + Complex_number cross_ratio; + for (int k=0; k> line; + index1 = std::stoi(line); + s >> line; + index2 = std::stoi(line); + dart_1 = darts_by_id[index1]; + dart_2 = darts_by_id[index2]; + _combinatorial_map.template sew<2>(dart_1, dart_2); + s >> cross_ratio; + _combinatorial_map.template set_attribute<1>(dart_1, _combinatorial_map.template create_attribute<1>(cross_ratio)); + } + + // (*) here + if (_has_anchor){ + _anchor.dart = darts_by_id[anchor_dart_id]; + } +} + +//////////////////////////////////////////////////////////////////////////////// + +template +std::ostream& operator<<(std::ostream& s, const Triangulation_on_hyperbolic_surface_2& triangulation){ + triangulation.to_stream(s); + return s; +} + +template +void operator>>(std::istream& s, Triangulation_on_hyperbolic_surface_2& triangulation){ + triangulation.from_stream(s); +} + +//////////////////////////////////////////////////////////////////////////////// + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_descriptor Triangulation_on_hyperbolic_surface_2::ccw(Dart_descriptor dart){ + return _combinatorial_map.beta(dart, 1); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_descriptor Triangulation_on_hyperbolic_surface_2::cw(Dart_descriptor dart){ + return _combinatorial_map.beta(dart, 0); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_descriptor Triangulation_on_hyperbolic_surface_2::opposite(Dart_descriptor dart){ + return _combinatorial_map.opposite(dart); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_const_descriptor Triangulation_on_hyperbolic_surface_2::const_ccw(Dart_const_descriptor dart) const{ + return _combinatorial_map.beta(dart, 1); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_const_descriptor Triangulation_on_hyperbolic_surface_2::const_cw(Dart_const_descriptor dart) const{ + return _combinatorial_map.beta(dart, 0); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Dart_const_descriptor Triangulation_on_hyperbolic_surface_2::const_opposite(Dart_const_descriptor dart) const{ + return _combinatorial_map.opposite(dart); +} + +//////////////////////////////////////////////////////////////////////////////// + +template +typename Triangulation_on_hyperbolic_surface_2::Complex_number Triangulation_on_hyperbolic_surface_2::get_cross_ratio(Dart_const_descriptor dart) const{ + return _combinatorial_map.template info_of_attribute<1>(_combinatorial_map.template attribute<1>(dart)); +} + +//////////////////////////////////////////////////////////////////////////////// + +template +void Triangulation_on_hyperbolic_surface_2::copy_from(Combinatorial_map_with_cross_ratios& cmap){ + //_combinatorial_map.copy_from_const(cmap); + _combinatorial_map.copy(cmap); + _has_anchor = false; +} + +template +void Triangulation_on_hyperbolic_surface_2::copy_from(Combinatorial_map_with_cross_ratios& cmap, const Anchor& anchor){ + // Because of the anchor, we must operate the copy ourself + _combinatorial_map.clear(); + + // Copy the triangles and fill the darts conversion table + std::map darts_table; + for (typename Face_const_range::const_iterator it=cmap.template one_dart_per_cell<2>().begin(); it!=cmap.template one_dart_per_cell<2>().end(); ++it){ + Dart_descriptor new_dart = _combinatorial_map.make_combinatorial_polygon(3); + darts_table[it] = new_dart; + darts_table[cmap.beta(it,0)] = _combinatorial_map.beta(new_dart,0); + darts_table[cmap.beta(it,1)] = _combinatorial_map.beta(new_dart,1); + } + + // Sew the edges and set their cross-ratios + for (typename Edge_const_range::const_iterator it=cmap.template one_dart_per_cell<1>().begin(); it!=cmap.template one_dart_per_cell<1>().end(); ++it){ + Dart_descriptor dart_1 = darts_table[it]; + Dart_descriptor dart_2 = darts_table[cmap.opposite(it)]; + Complex_number cratio = cmap.template info_of_attribute<1>(cmap.template attribute<1>(it)); + + _combinatorial_map.template sew<2>(dart_1, dart_2); + _combinatorial_map.template set_attribute<1>(dart_1, _combinatorial_map.template create_attribute<1>(cratio)); + } + + cmap.opposite(anchor.dart); + + // Set the anchor + _anchor.dart = darts_table[anchor.dart]; + for (int k=0; k<3; k++){ + _anchor.vertices[k] = anchor.vertices[k]; + } + _has_anchor = true; +} + +//////////////////////////////////////////////////////////////////////////////// + +template +typename Traits::Complex Triangulation_on_hyperbolic_surface_2::cross_ratio(const Point& a, const Point& b, const Point& c, const Point& d) const{ + Complex_number za (a.x(), a.y()); + Complex_number zb (b.x(), b.y()); + Complex_number zc (c.x(), c.y()); + Complex_number zd (d.x(), d.y()); + return (zd-zb)*(zc-za) / ((zd-za)*(zc-zb)); +} + +template +typename Triangulation_on_hyperbolic_surface_2::Point Triangulation_on_hyperbolic_surface_2::fourth_point_from_cross_ratio(const Point& a, const Point& b, const Point& c, const Complex_number& cratio) const{ + Complex_number za (a.x(), a.y()); + Complex_number zb (b.x(), b.y()); + Complex_number zc (c.x(), c.y()); + Complex_number result = ( cratio*za*(zc-zb) + zb*(za-zc) ) / ( cratio*(zc-zb) + (za-zc) ); + return Point(result.real(), result.imag()); +} + +} // namespace CGAL + +#endif // CGAL_TRIANGULATION_ON_HYPERBOLIC_SURFACE_2