The traits class now uses circular kernel.

- dual functions in Delaunay class fixed (but code should still be improved)
- examples and demos modified accordingly
- broken: Hyperbolic_random_points_in_disc_2 and benchmarks (problems with number types)
This commit is contained in:
Monique Teillaud 2016-08-11 17:15:18 +02:00 committed by Iordan Iordanov
parent 8c88155cc7
commit bd9b60f7cc
11 changed files with 420 additions and 343 deletions

View File

@ -1,28 +1,55 @@
==== new ====
add operator == for triangulations
add operator == for triangulations (?)
** clean approximations ---> use CK
** traits bisectors
the construction of the (hyperbolic) segment should not be approximated in the traits
the approximation should be done for the demo only
class Construct_hyperbolic_bisector_2
Hyperbolic_segment_2 operator()(Point_2 p, Point_2 q, Point_2 r)
at the end
the following lines are giving the wrong arc (wrong orientation), aren't they?
CGAL_triangulation_assertion(assign(pair,inters[1]));
if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE )
return Circular_arc_2( *c_pq, pair.first, a);
return Circular_arc_2( *c_pq, a, pair.first);
but this is what works on the demo...
understand why
similar: clean find_intersection
(using the circular kernel will do)
remove variant for supporting circle or line of bisector
call it only when we know that it is a circle
it will simplyfy the code of Construct_hyperbolic_bisector_2 at least in some cases
fix construct_circumcenter. no approx allowed
(again, CK will do)
test bisectors dual functions in special cases of euclidean line segments
CK, or at least Cartesian with sqrt (?)
** Hyperbolic_random_points_in_disc
repair
number types have changed (using CK)
** benchmarks
idem, repair number types
** doc
** tests
** test suite
** only keep relevant files for the submission
========== demo
--- input from file does not do anything...
--- clean files
do we need demo/Hyperbolic_triangulation_2/include/* for the basic demo?
todo diff with similar .h in Graphics view (after updating the branch)
--- clean buttons
remove translations, add circumcircle, etc
and put appropriate window title
--- fix aretes delaunay presque rectilineaires pas au bon endroit
--- create an adapted Graphicsitem that does not require finite_* stuff. Replace by hyperbolic_*
=> understand apply_to_range
the graphicsitem should also not require Segment_2 or Line_segment_2 (does it?)
@ -60,13 +87,6 @@ CGAL::https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/Unique_sqrt_ex
Done. std::sqrt is replaced by CGAL::sqrt everywhere,
because it's used only for computation of the Voronoi diagram.
--- use the Circular_kernel?
because Do_intersect is in the Circular_kernel
Done. CGAL::do_intersect is a solution.
No need for the Circular_kernel.
--> (MT) wrong for constrcutions. Approximations are used in several places.
========== future
--- make Constrained_delaunay_... work with Hyperbolic_triangulation_2

View File

@ -3,12 +3,12 @@
#include <CGAL/Hyperbolic_random_points_in_disc_2.h>
// CGAL headers
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
#include <CGAL/Timer.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_circular_kernel_2 K;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2<K> Gt;
typedef K::Point_2 Point_2;

View File

@ -3,12 +3,12 @@
#include <CGAL/Hyperbolic_random_points_in_disc_2.h>
// CGAL headers
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
#include <CGAL/Timer.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_circular_kernel_2 K;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2<K> Gt;
typedef K::Point_2 Point_2;

View File

@ -1,14 +1,12 @@
#include <fstream>
// CGAL headers
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
// to be deleted
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Qt/HyperbolicPainterOstream.h>
//
#include <CGAL/point_generators_2.h>
@ -22,13 +20,13 @@
// GraphicsView items and event filters (input classes)
#include "TriangulationCircumcircle.h"
#include "TriangulationMovingPoint.h"
#include "TriangulationConflictZone.h"
#include "TriangulationRemoveVertex.h"
#include "TriangulationPointInputAndConflictZone.h"
#include <CGAL/Qt/TriangulationGraphicsItem.h>
#include <CGAL/Qt/VoronoiGraphicsItem.h>
#include <CGAL/Qt/HyperbolicVoronoiGraphicsItem.h>
// for viewportsBbox
#include <CGAL/Qt/utility.h>
@ -37,7 +35,8 @@
#include "ui_Hyperbolic_translations_2.h"
#include <CGAL/Qt/DemosMainWindow.h>
typedef CGAL::Exact_predicates_exact_constructions_kernel R;
//typedef CGAL::Exact_predicates_exact_constructions_kernel R;
typedef CGAL::Exact_circular_kernel_2 R;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2<R> K;
typedef K::Point_2 Point_2;

View File

@ -1,15 +1,12 @@
#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/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
// to be deleted
#include <CGAL/Qt/HyperbolicPainterOstream.h>
//
#include <CGAL/point_generators_2.h>
@ -31,7 +28,8 @@
#include "TriangulationRemoveVertex.h"
#include "TriangulationPointInputAndConflictZone.h"
//#include <CGAL/Qt/TriangulationGraphicsItem.h>
#include <CGAL/Qt/VoronoiGraphicsItem.h>
#include <CGAL/Qt/HyperbolicVoronoiGraphicsItem.h>
// store color
#include <TranslationInfo.h>
@ -55,8 +53,7 @@
#include <CGAL/Qt/DemosMainWindow.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel InR;
typedef CGAL::Exact_predicates_exact_constructions_kernel R;
typedef CGAL::Exact_circular_kernel_2 R;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2<R> K;
typedef K::Point_2 Point_2;

View File

@ -3,7 +3,7 @@
// CGAL headers
#include <CGAL/IO/io.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_triangulation_face_base_with_info_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
@ -11,7 +11,7 @@
#include <CGAL/IO/Color.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_circular_kernel_2 K;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2< K > Gt;
typedef Gt::Point_2 Point_2;

View File

@ -3,8 +3,7 @@
// CGAL headers
#include <CGAL/IO/io.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_2.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
@ -13,7 +12,7 @@
#include <CGAL/point_generators_2.h>
#include <CGAL/Timer.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_circular_kernel_2 K;
typedef CGAL::Hyperbolic_Delaunay_triangulation_traits_2<K> Gt;
typedef Gt::Point_2 Point_2;
@ -25,31 +24,31 @@ typedef CGAL::Hyperbolic_Delaunay_triangulation_2<Gt> Dt;
int main()
{
CGAL::Timer timer;
// typedef CGAL::Creator_uniform_2<FT, Point_2> Creator;
typedef CGAL::Creator_uniform_2<FT, Point_2> Creator;
// CGAL::Random_points_in_disc_2<Point_2, Creator> in_disc(FT(1));
CGAL::Random_points_in_disc_2<Point_2, Creator> in_disc(FT(1));
// int n = 10;
// std::cout << "Number of points: " << n << std::endl;
int n = 100;
std::cout << "Number of points: " << n << std::endl;
// std::vector<Point_2> pts(n);
// std::vector<Point_2>::iterator ip;
// // Generating n random points
// for (int i=0 ; i < n ; i++) {
// pts.at(i) = *in_disc;
// in_disc++;
// }
std::vector<Point_2> pts;
std::vector<Point_2> pts(n);
std::vector<Point_2>::iterator ip;
Point_2 p;
std::ifstream ifs("input-file");
while(ifs >> p) {
pts.push_back(p);
// Generating n random points
for (int i=0 ; i < n ; i++) {
pts.at(i) = *in_disc;
in_disc++;
}
std::cout << "number of points " << std::distance(pts.begin(),pts.end()) << std::endl << std::endl;
// std::vector<Point_2> pts;
// std::vector<Point_2>::iterator ip;
// Point_2 p;
// std::ifstream ifs("input-file");
// while(ifs >> p) {
// pts.push_back(p);
// }
// std::cout << "number of points " << std::distance(pts.begin(),pts.end()) << std::endl << std::endl;
std::cout << "check for hyperbolic faces during insertion" << std::endl;

View File

@ -64,8 +64,8 @@ public:
typedef Gt Geom_traits;
typedef typename Geom_traits::FT FT;
typedef typename Geom_traits::Point_2 Point;
typedef typename Geom_traits::Hyperbolic_segment_2 Segment;
typedef typename Geom_traits::Circular_arc_point_2 Circular_arc_point;
typedef typename Geom_traits::Hyperbolic_segment_2 Hyperbolic_segment;
Hyperbolic_Delaunay_triangulation_2(const Gt& gt = Gt())
: Delaunay_triangulation_2<Gt,Tds>(gt) {}
@ -477,7 +477,7 @@ public:
Finite_edges_iterator finite_edges_begin() const { return hyperbolic_edges_begin(); }
Finite_edges_iterator finite_edges_end() const { return hyperbolic_edges_end(); }
Point
Circular_arc_point
dual(Face_handle f) const
{
CGAL_triangulation_precondition (!this->is_non_hyperbolic(f));
@ -486,65 +486,70 @@ public:
( f->vertex(0)->point(), f->vertex(1)->point(), f->vertex(2)->point());
}
Object
Hyperbolic_segment
dual(const Edge& e) const
{
return dual(e.first, e.second);
}
Object
Hyperbolic_segment
dual(Face_handle f, int i) const
{
CGAL_triangulation_precondition (!this->is_non_hyperbolic(f,i));
if(this->dimension() == 1) {
const Point& p = f->vertex(cw(i))->point();
const Point& q = f->vertex(ccw(i))->point();
Point p = f->vertex(cw(i))->point();
Point q = f->vertex(ccw(i))->point();
// hyperbolic line
Segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q);
return make_object(line);
Hyperbolic_segment line = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q);
return line;
}
Face_handle n = f->neighbor(i);
int in = n->index(f);
//TODO MT store values of bools to avoid recomputing is-hyperbolic several times
// boths faces are non_hyperbolic, but the incident edge is hyperbolic
if(is_non_hyperbolic(f) && is_non_hyperbolic(n)){
const Point& p = f->vertex(cw(i))->point();
const Point& q = f->vertex(ccw(i))->point();
if( is_non_hyperbolic(f) && is_non_hyperbolic(n) ){
const Point& p = f->vertex(ccw(i))->point();
const Point& q = f->vertex(cw(i))->point();
// hyperbolic line
Segment line =
Hyperbolic_segment line =
this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q);
return make_object(line);
return line;
}
// both faces are finite
if(!is_non_hyperbolic(f) && !is_non_hyperbolic(n)) {
Segment s =
this->geom_traits().construct_segment_2_object()(dual(f),dual(n));
return make_object(s);
// both faces are hyperbolic
if( !is_non_hyperbolic(f) && !is_non_hyperbolic(n) ) {
const Point& p = f->vertex(ccw(i))->point();
const Point& q = f->vertex(cw(i))->point();
Hyperbolic_segment s =
this->geom_traits().construct_hyperbolic_bisector_2_object()
(p,q,f->vertex(i)->point(),n->vertex(in)->point());
//TODO MT cut edge at dual points !!!!
return s;
}
// one of the incident faces is non_hyperbolic
Face_handle finite_face = f;
Face_handle hyp_face = f;
if(is_non_hyperbolic(f)) {
finite_face = n;
hyp_face = n;
i = in;
}
const Point& p = finite_face->vertex(cw(i))->point();
const Point& q = finite_face->vertex(ccw(i))->point();
const Point& p = hyp_face->vertex(ccw(i))->point();
const Point& q = hyp_face->vertex(cw(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);
Hyperbolic_segment ray = this->geom_traits().construct_hyperbolic_bisector_2_object()(p,q,hyp_face->vertex(i)->point());
// TODO MT cut edge at dual point !!!
// Segment ray = this->geom_traits().construct_ray_2_object()(dual(finite_face), line);
return ray;
}
};

View File

@ -22,7 +22,7 @@
#ifndef CGAL_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H
#define CGAL_HYPERBOLIC_DELAUNAY_TRIANGULATION_TRAITS_2_H
#include <CGAL/Regular_triangulation_filtered_traits_2.h>
#include <CGAL/Regular_triangulation_euclidean_traits_2.h>
#include "boost/tuple/tuple.hpp"
#include "boost/variant.hpp"
@ -31,19 +31,25 @@ namespace CGAL {
template < class R >
class Hyperbolic_Delaunay_triangulation_traits_2
: public R
// R is supposed to be a model of CircularKernel2
{
public:
typedef typename R::FT FT;
typedef typename R::Point_2 Point_2;
typedef typename R::Circle_2 Circle_2;
typedef typename R::Line_2 Euclidean_line_2;
typedef boost::variant<Circle_2,Euclidean_line_2> Euclidean_circle_or_line_2;
typedef boost::tuple<Circle_2, Point_2, Point_2> Arc_2;
typedef typename R::Circular_arc_2 Circular_arc_2;
typedef typename R::Line_arc_2 Line_arc_2;
typedef typename R::Circular_arc_point_2 Circular_arc_point_2;
typedef typename R::Segment_2 Euclidean_segment_2; //only used internally here
typedef boost::variant<Arc_2, Euclidean_segment_2> Hyperbolic_segment_2;
typedef boost::variant<Circular_arc_2, Line_arc_2> Hyperbolic_segment_2;
typedef typename R::Compare_x_2 Compare_x_2;
typedef typename R::Compare_y_2 Compare_y_2;
typedef typename R::Compare_distance_2 Compare_distance_2;
typedef typename R::Orientation_2 Orientation_2;
typedef typename R::Side_of_oriented_circle_2 Side_of_oriented_circle_2;
@ -57,12 +63,7 @@ public:
typedef typename R::Construct_bisector_2 Construct_Euclidean_bisector_2;
typedef typename R::Construct_midpoint_2 Construct_Euclidean_midpoint_2;
typedef typename R::Compute_squared_distance_2 Compute_squared_Euclidean_distance_2;
typedef typename R::Line_2 Euclidean_line_2;
typedef typename R::Vector_2 Vector_2;
// used by Is_hyperbolic
typedef typename R::Vector_3 Vector_3;
typedef typename R::Point_3 Point_3;
typedef typename R::Has_on_bounded_side_2 Has_on_bounded_side_2;
// MT useless?
// typedef Hyperbolic_Delaunay_triangulation_traits_2<R> Self;
// typedef typename R::RT RT;
@ -74,9 +75,8 @@ public:
// typedef typename R::Iso_rectangle_2 Iso_rectangle_2;
// typedef typename R::Angle_2 Angle_2;
// typedef typename R::Less_x_2 Less_x_2;
// typedef typename R::Less_y_2 Less_y_2;
// typedef typename R::Compare_distance_2 Compare_distance_2;
typedef typename R::Less_x_2 Less_x_2;
typedef typename R::Less_y_2 Less_y_2;
// typedef typename R::Construct_triangle_2 Construct_triangle_2;
// typedef typename R::Construct_direction_2 Construct_direction_2;
@ -84,15 +84,14 @@ public:
class Construct_hyperbolic_segment_2
{
typedef typename CGAL::Regular_triangulation_filtered_traits_2<R> Regular_geometric_traits_2;
typedef typename CGAL::Regular_triangulation_euclidean_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_hyperbolic_segment_2()
{
}
Construct_hyperbolic_segment_2()
{}
Hyperbolic_segment_2 operator()(const Point_2& p, const Point_2& q) const
{
@ -106,19 +105,19 @@ public:
Weighted_point_2 wo(Point_2(o), FT(1)); // Poincaré circle
Bare_point center = Construct_weighted_circumcenter_2()(wp, wo, wq);
FT radius = Compute_squared_Euclidean_distance_2()(p, center);
FT sq_radius = Compute_squared_Euclidean_distance_2()(p, center);
Circle_2 circle( center, radius);
Circle_2 circle(center, sq_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 Circular_arc_2(circle, p, q);
}
return Arc_2(circle, q, p);
return Circular_arc_2(circle, q, p);
}
};
}; // end Construct_hyperbolic_segment_2
Construct_hyperbolic_segment_2
construct_hyperbolic_segment_2_object() const
@ -134,38 +133,71 @@ public:
{
public:
// TODO: improve this function
Point_2 operator()(Point_2 p, Point_2 q, Point_2 r)
{
CGAL_triangulation_assertion_code(Origin oo; Point_2 poo(oo); Circle_2 co(poo,FT(1)));
CGAL_triangulation_assertion(co.bounded_side(p) == ON_BOUNDED_SIDE);
CGAL_triangulation_assertion(co.bounded_side(q) == ON_BOUNDED_SIDE);
CGAL_triangulation_assertion(co.bounded_side(r) == ON_BOUNDED_SIDE);
Circle_2 circle(p, q, r);
// circle must be inside the unit one
CGAL_triangulation_assertion(do_intersect(co, circle) == false);
Circular_arc_point_2 operator()(Point_2 p, Point_2 q, Point_2 r)
{
Origin o;
Point_2 po = Point_2(o);
if (circle.center() == po)
{ return po; }
Circle_2 l_inf(po, FT(1));
Euclidean_circle_or_line_2 bis_pq = Construct_circle_or_line_supporting_bisector()(p,q);
Euclidean_circle_or_line_2 bis_qr = Construct_circle_or_line_supporting_bisector()(q,r);
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() + FT(1); // Poincare disc has radius 1
FT D = b*b - 4*a;
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;
if ( Compare_distance_2()(po,p,q) == EQUAL &&
Compare_distance_2()(po,p,r) == EQUAL )
return po;
// now supporting objects cannot both be Euclidean lines
std::pair<Circular_arc_point_2, unsigned> pair;
Euclidean_line_2* l;
Circle_2* c;
if ( Circle_2* c_pq = boost::get<Circle_2>(&bis_pq) )
{
if ( Circle_2* c_qr = boost::get<Circle_2>(&bis_qr) )
{
typedef typename CK2_Intersection_traits<R, Circle_2, Circle_2>::type Intersection_result;
std::vector< Intersection_result > inters;
intersection(*c_pq, *c_qr, std::back_inserter(inters));
CGAL_triangulation_assertion(assign(pair,inters[0]));
if ( pair.second == 1 )
{
if ( Has_on_bounded_side_2()( l_inf, pair.first ) )
return pair.first;
CGAL_triangulation_assertion(assign(pair,inters[1]));
return pair.first;
}
return pair.first;
}
// here bis_qr is a line
l = boost::get<Euclidean_line_2>(&bis_qr);
c = c_pq;
}
// here bis_pq is a line
l = boost::get<Euclidean_line_2>(&bis_pq);
c = boost::get<Circle_2>(&bis_qr);
typedef typename CK2_Intersection_traits<R, Euclidean_line_2, Circle_2>::type Intersection_result;
std::vector< Intersection_result > inters;
intersection(*l, *c, std::back_inserter(inters));
CGAL_triangulation_assertion(assign(pair,inters[0]));
if ( pair.second == 1 )
{
if ( Has_on_bounded_side_2()( l_inf, pair.first ) )
return pair.first;
CGAL_triangulation_assertion(assign(pair,inters[1]));
return pair.first;
}
return pair.first;
}
};
}; // end Construct_hyperbolic_circumcenter_2
Hyperbolic_Delaunay_triangulation_traits_2()
{}
@ -202,120 +234,141 @@ public:
class Construct_hyperbolic_bisector_2
{
public:
Construct_hyperbolic_bisector_2()
{}
Construct_hyperbolic_bisector_2()
{}
// constructs a hyperbolic line
Hyperbolic_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_Euclidean_distance_2 dist = Compute_squared_Euclidean_distance_2();
Origin o;
Point_2 po = Point_2(o);
FT dif = dist(po, p) - dist(po, q);
FT eps = 0.0000000001;
Circle_2 l_inf = Circle_2(po,FT(1));
// Bisector is straight in euclidean sense
if(dif > -eps && dif < eps){
// ideally
//if(Compare_distance_2()(origin, p, q) == EQUAL){
// TODO: calling R::Construct_bisector
Euclidean_line_2 l = Construct_Euclidean_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 Euclidean_segment_2(points.first, points.second);
}
return Euclidean_segment_2(points.second, points.first);
if ( Compare_distance_2()(po,p,q) == EQUAL ){
Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p,q);
if ( Less_y_2()(p,q) )
{ return Line_arc_2( l, l_inf, false, l_inf, true ); }
return Line_arc_2( l, l_inf, true, l_inf, false );
}
Circle_2 c = construct_supporting_circle_of_bisector(p, q);
// compute the ending points
std::pair<Point_2, Point_2> points = find_intersection(c);
Euclidean_circle_or_line_2
bis_pq = Construct_circle_or_line_supporting_bisector()(p, q);
Circle_2* c = boost::get<Circle_2>(&bis_pq);
if(Orientation_2()(points.first, points.second, c.center()) == LEFT_TURN) {
return Arc_2(c, points.first, points.second);
if ( Less_y_2()(po,c->center()) )
{ return Circular_arc_2( *c, l_inf, true, l_inf, false ); }
else if ( Less_y_2()(c->center(),po) )
{ return Circular_arc_2( *c, l_inf, false, l_inf, true ); }
// the center of the circle is on the x-axis
if ( Less_x_2()(po,c->center()) )
{ return Circular_arc_2( *c, l_inf, true, l_inf, false ); }
return Circular_arc_2( *c, l_inf, false, l_inf, true);
}
// constructs the hyperbolic bisector of segment [p,q] limited by
// circumcenter(p,q,r) on one side
// and circumcenter(p,s,q) on the other side
Hyperbolic_segment_2
operator()(Point_2 p, Point_2 q, Point_2 r, Point_2 s)
{
CGAL_triangulation_precondition
( (Orientation_2()(p,q,r) == ON_POSITIVE_SIDE)
&& (Orientation_2()(p,s,q) == ON_POSITIVE_SIDE) );
CGAL_triangulation_precondition
( (Side_of_oriented_circle_2()(p,q,r,s) == ON_NEGATIVE_SIDE)
&& (Side_of_oriented_circle_2()(p,s,q,r) == ON_NEGATIVE_SIDE) );
Origin o;
Point_2 po = Point_2(o);
// TODO MT this is non-optimal...
// the bisector is already computed here
// and it will be recomputed below
Circular_arc_point_2 a = Construct_hyperbolic_circumcenter_2()(p,q,r);
Circular_arc_point_2 b = Construct_hyperbolic_circumcenter_2()(p,s,q);
if ( Compare_distance_2()(po, p, q) == EQUAL ){
Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p,q);
return Line_arc_2(l,a,b);
}
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_of_bisector(Point_2 p, Point_2 q) const
Euclidean_circle_or_line_2
bis_pq = Construct_circle_or_line_supporting_bisector()(p, q);
Circle_2* c_pq = boost::get<Circle_2>(&bis_pq);
if ( Compare_distance_2()(po,p,q) == POSITIVE )
// then p is inside the supporting circle
{ return Circular_arc_2(*c_pq,b,a);}
return Circular_arc_2(*c_pq,a,b);
}
// constructs the hyperbolic bisector of segment [p,q]
// limited by hyperbolic circumcenter(p,q,r) on one side
// and going to the infinite line on the other side
Hyperbolic_segment_2 operator()(Point_2 p, Point_2 q, Point_2 r)
{
// 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 = (FT(1) - oq) / (op - oq); // Poincare disc has radius 1
FT x = alpha*p.x() + (1-alpha)*q.x();
FT y = alpha*p.y() + (1-alpha)*q.y();
FT radius = x*x + y*y - FT(1);
//improve
Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p, q);
Point_2 middle = Construct_Euclidean_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);
CGAL_triangulation_precondition
( Orientation_2()(p,q,r) == POSITIVE );
Origin o;
Point_2 po = Point_2(o);
Circle_2 l_inf(po, FT(1));
// TODO MT this is non-optimal...
// the bisector is computed (at least) twice
Circular_arc_point_2 a = Construct_hyperbolic_circumcenter_2()(p,q,r);
if ( Compare_distance_2()(po, p, q) == EQUAL ){
Euclidean_line_2 bis_pq = Construct_Euclidean_bisector_2()(p,q);
typedef typename
CK2_Intersection_traits<R, Euclidean_line_2, Circle_2>::type
Intersection_result;
std::vector< Intersection_result > inters;
intersection(bis_pq, l_inf, std::back_inserter(inters));
std::pair<Circular_arc_point_2, unsigned> pair;
CGAL_triangulation_assertion(assign(pair,inters[0]));
CGAL_triangulation_assertion(pair.second == 1);
if ( Less_y_2()(p,q) )
return Line_arc_2(bis_pq,a,pair.first);
CGAL_triangulation_assertion(assign(pair,inters[1]));
CGAL_triangulation_assertion(pair.second == 1);
return Line_arc_2(bis_pq,a,pair.first);
}
Euclidean_circle_or_line_2
bis_pq = Construct_circle_or_line_supporting_bisector()(p, q);
Circle_2* c_pq = boost::get<Circle_2>(&bis_pq);
return Circle_2(Point_2(x, y), radius, COUNTERCLOCKWISE);
Point_2 approx_a(to_double(a.x()),to_double(a.y()));
typedef typename
CK2_Intersection_traits<R, Circle_2, Circle_2>::type
Intersection_result;
std::vector< Intersection_result > inters;
intersection(*c_pq, l_inf, std::back_inserter(inters));
std::pair<Circular_arc_point_2, unsigned> pair;
CGAL_triangulation_assertion(assign(pair,inters[0]));
CGAL_triangulation_assertion(pair.second == 1);
Point_2 approx_pinf(to_double(pair.first.x()), to_double(pair.first.y()));
Point_2 approx_c(to_double(c_pq->center().x()),
to_double(c_pq->center().y()));
if ( Orientation_2()(p,q,approx_pinf) == NEGATIVE ) {
if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE )
return Circular_arc_2( *c_pq, a, pair.first );
return Circular_arc_2( *c_pq, pair.first, a);
}
CGAL_triangulation_assertion(assign(pair,inters[1]));
if ( Orientation_2()(approx_c,approx_a,approx_pinf) == POSITIVE )
return Circular_arc_2( *c_pq, pair.first, a);
return Circular_arc_2( *c_pq, a, pair.first);
}
// Find intersection of an input circle orthogonal to the Poincare 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; */
FT b = -x;
FT c = 1-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 = (FT(1) - x1*x)/y;
FT y2 = (FT(1) - x2*x)/y;
return std::make_pair(Point_2(x1, y1), Point_2(x2, y2));
}
// Find intersection of an input line orthogonal to the Poincare 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 = FT(1)/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);
}
};
}; // end Construct_hyperbolic_bisector_2
Construct_hyperbolic_bisector_2
construct_hyperbolic_bisector_2_object() const
@ -325,45 +378,13 @@ public:
construct_Euclidean_bisector_2_object() const
{ return Construct_Euclidean_bisector_2(); }
class Construct_ray_2
{
public:
Construct_ray_2()
{}
Hyperbolic_segment_2 operator()(Point_2 p, Hyperbolic_segment_2 l) const
{
if(Euclidean_segment_2* s = boost::get<Euclidean_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 Hyperbolic_segment_2();
}
Hyperbolic_segment_2 operator()(Point_2 p, Euclidean_segment_2 s) const
{
return Euclidean_segment_2(p, s.target());
}
};
Construct_ray_2
construct_ray_2_object() const
{ return Construct_ray_2(); }
// For details see the JoCG paper (5:56-85, 2014)
class Is_hyperbolic
{
public:
typedef typename R::Vector_3 Vector_3;
typedef typename R::Point_3 Point_3;
bool operator() (const Point_2& p0, const Point_2& p1, const Point_2& p2) const
{
Vector_3 v0 = Vector_3(p0.x()*p0.x() + p0.y()*p0.y(),
@ -425,13 +446,71 @@ public:
return 1;
}
};
}; // end Is_hyperbolic
Is_hyperbolic
Is_hyperbolic_object() const
{ return Is_hyperbolic(); }
// do not document
// constructs the Euclidean circle or line supporting the hyperbolic
// bisector of two points
class Construct_circle_or_line_supporting_bisector
{
public:
Construct_circle_or_line_supporting_bisector()
{}
Euclidean_circle_or_line_2
operator()(Point_2 p, Point_2 q) const
{
Origin o;
Point_2 po = Point_2(o);
typedef typename R::Point_3 Point_3;
if ( Compare_distance_2()(po,p,q) == EQUAL )
{ return Construct_Euclidean_bisector_2()(p,q); }
FT dop2 = p.x()*p.x() + p.y()*p.y();
FT doq2 = q.x()*q.x() + q.y()*q.y();
Point_3 p3( p.x(), p.y(), dop2 );
Point_3 q3( q.x(), q.y(), doq2 );
// TODO MT improve
// The cirle belongs to the pencil with limit points p and q
// 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 = 1 (= radius of the Poincare disc)
FT op = p.x()*p.x() + p.y()*p.y();
FT oq = q.x()*q.x() + q.y()*q.y();
FT alpha = (FT(1) - oq) / (op - oq);
FT x = alpha*p.x() + (1-alpha)*q.x();
FT y = alpha*p.y() + (1-alpha)*q.y();
FT sq_radius = x*x + y*y - FT(1);
// TODO MT improve
// ?? orientation should depend on
// Compare_distance(O,p,q)
// so that p always on positive side
// ???
// CK does not care about orientation, circular arcs are
// considered in CCW order in any case
Euclidean_line_2 l = Construct_Euclidean_bisector_2()(p, q);
Point_2 middle = Construct_Euclidean_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), sq_radius, CLOCKWISE); }
return Circle_2(Point_2(x, y), sq_radius, COUNTERCLOCKWISE);
}
}; // end Construct_supporting_circle_of_bisector
};
// Take out the code below to some separate file
#ifdef CGAL_EXACT_PREDICATES_EXACT_CONSTRUCTIONS_KERNEL_H

View File

@ -16,79 +16,65 @@
// $Id:
//
//
// Author(s) : Mikhail Bogdanov
// Author(s) : Monique Teillaud <Monique.Teillaud@inria.fr>
// Mikhail Bogdanov
#ifndef CGAL_HYPERBOLIC_PAINTER_OSTREAM_H
#define CGAL_HYPERBOLIC_PAINTER_OSTREAM_H
#include <CGAL/Qt/PainterOstream.h>
#include <CGAL/Hyperbolic_Delaunay_triangulation_traits_2.h>
namespace CGAL{
namespace Qt {
template <typename K>
class PainterOstream<Hyperbolic_Delaunay_triangulation_traits_2<K> > : public PainterOstream<K> {
class PainterOstream<Hyperbolic_Delaunay_triangulation_traits_2<K> >
: public PainterOstream<K>
{
public:
typedef PainterOstream<K> Base;
typedef PainterOstream<Hyperbolic_Delaunay_triangulation_traits_2<K> > Self;
typedef Hyperbolic_Delaunay_triangulation_traits_2<K> Gt;
typedef PainterOstream<Gt> Self;
typedef typename Gt::Hyperbolic_segment_2 Hyperbolic_segment_2;
typedef typename Gt::Point_2 Point_2;
typedef Line_arc_2<K> Line_arc_2;
typedef Circular_arc_2<K> Circular_arc_2;
typedef Circular_arc_point_2<K> Circular_arc_point_2;
private:
typedef typename Gt::Segment_2 Segment_2;
typedef typename Gt::Line_segment_2 Line_segment_2;
typedef typename Gt::Arc_2 Arc_2;
//typedef typename Gt::Line_2 Line_2;
typedef typename K::Point_2 Point_2;
typedef typename K::Circle_2 Circle_2;
public:
PainterOstream(QPainter* p, QRectF rect = QRectF())
: Base(p, rect), qp(p), convert(rect) {}
: Base(p, rect), qp(p), convert(rect)
{}
using Base::operator <<;
PainterOstream& operator << (const Segment_2& s)
{
if(const Line_segment_2* line_segment = boost::get<Line_segment_2>(&s)){
operator << (*line_segment);
return *this;
PainterOstream& operator << (Hyperbolic_segment_2 s)
{
if (const Line_arc_2* seg = boost::get<Line_arc_2>(&s)) {
operator << (*seg);
return *this;
}
Circular_arc_2* arc = boost::get<Circular_arc_2>(&s);
if (arc->squared_radius() > 10 )
// due to rounding, the arc drawn does not look like it
// passes through the endpoints
// so we replace the arc by a line segment
{
Point_2 source(to_double(arc->source().x()),to_double(arc->source().y()));
Point_2 target(to_double(arc->target().x()),to_double(arc->target().y()));
const Line_arc_2 seg(source,target);
operator << (seg);
return *this;
}
operator << (*arc);
return *this;
}
if(const Arc_2* arc = boost::get<Arc_2>(&s)){
const Circle_2& circle = get<0>(*arc);
const Point_2& center = circle.center();
const Point_2& source = get<1>(*arc);
const Point_2& target = get<2>(*arc);
if (circle.squared_radius() > 10) {
const Line_segment_2 seg(source,target);
operator << (seg);
return *this;
}
double asource = std::atan2( -to_double(source.y() - center.y()),
to_double(source.x() - center.x()));
double atarget = std::atan2( -to_double(target.y() - center.y()),
to_double(target.x() - center.x()));
std::swap(asource, atarget);
double aspan = atarget - asource;
if(aspan < 0.)
aspan += 2 * CGAL_PI;
const double coeff = 180*16/CGAL_PI;
this->qp->drawArc(this->convert(circle.bbox()),
(int)(asource * coeff),
(int)(aspan * coeff));
}
return *this;
}
private:
// ToDo: These objects must be deleted

View File

@ -108,20 +108,12 @@ VoronoiGraphicsItem<DT>::paint(QPainter *painter, const QStyleOptionGraphicsItem
for(typename DT::Finite_edges_iterator eit = dt->finite_edges_begin();
eit != dt->finite_edges_end();
eit++){
CGAL::Object o = dt->dual(*eit);
typename DT::Segment s;
typename DT::Geom_traits::Ray_2 r;
typename DT::Geom_traits::Line_2 l;
if(CGAL::assign(s,o)){
eit++)
{
typename DT::Hyperbolic_segment s = dt->dual(*eit);
pos << s;
} else if(CGAL::assign(r,o)) {
pos << r;
}else if(CGAL::assign(l,o)) {
pos << l;
}
}
// delete
painter->setPen(old);
}