diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/intersect.spheres b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/intersect.spheres index e6712c6ec6c..1d49fc29aea 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/intersect.spheres +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/intersect.spheres @@ -1,2 +1,2 @@ -0 0 0 1.2 1 -.2 1 1 1 1 +0 0 0 1 1 +0 1 1 1 1 diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/interactive_ray_shooting.cpp b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/interactive_ray_shooting.cpp index a6543e924d4..da351792f9b 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/interactive_ray_shooting.cpp +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/interactive_ray_shooting.cpp @@ -1,7 +1,7 @@ #define CGAL_CHECK_EXPENSIVE #define CGAL_CHECK_EXACTNESS -#include +#include #include #include #include @@ -157,7 +157,7 @@ int main(int argc, char *argv[]){ Do_work dw(atof(argv[1]), argv[2]); - Qt_debug_viewer_2 qtd(dw, argc, argv); + Qt_multithreaded_examiner_viewer_2 qtd(dw, argc, argv); return qtd(); } diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/makefile b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/makefile index b799590841a..c05d5c14120 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/makefile +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/makefile @@ -37,7 +37,9 @@ LDFLAGS = \ # target entries #---------------------------------------------------------------------# -all: mylibs interactive_point_insertion slice_kds show_slice show_circles interactive_ray_shooting interactive_point_location +all: mylibs interactive_point_insertion slice_kds show_slice show_circles interactive_ray_shooting + +#interactive_point_location #interactive_ray_shooting interactive_point_location random_ray_shoot random_point_locate point_locate diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/show_circles.cpp b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/show_circles.cpp index cb05e9fc17f..ed6ff4b4e52 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/show_circles.cpp +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/show_circles.cpp @@ -1,7 +1,7 @@ -#define CGAL_CHECK_EXPENSIVE -#define CGAL_CHECK_EXACTNESS +//#define CGAL_CHECK_EXPENSIVE +//#define CGAL_CHECK_EXACTNESS -#include +#include #include #include #include @@ -82,7 +82,7 @@ int main(int argc, char *argv[]){ } - app.setMainWidget( qtd); + app.setMainWidget( qtd->window()); qtd->show_everything(); qtd->show(); return app.exec(); diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice.h index cdc4fa12e70..e00e3244ec1 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice.h @@ -206,6 +206,10 @@ struct Slice { void audit_edge_collapse_event(Halfedge_handle h, Event_key ek); + void process_equal_centers_two_sphere_event(T::Key k, T::Key l, bool f); + void process_aligned_centers_two_sphere_event(T::Key k, T::Key l, bool f, + CGAL::Comparison_result c[2]); + // Functions to update certificates ---------------------------------- void clean_edge(Halfedge_handle h); @@ -342,10 +346,10 @@ struct Slice { */ - Face_handle locate_point(const T::Sphere_point_3 &ep, T::Key ind); + Face_handle locate_point(const T::Sphere_point_3 &ep); template - Face_handle locate_point(It b, It e, const T::Sphere_point_3 &ep, T::Key ind); + Face_handle locate_point(It b, It e, const T::Sphere_point_3 &ep); Halfedge_handle shoot_rule(const T::Sphere_point_3& source, Face_handle f, @@ -364,29 +368,24 @@ struct Slice { //enum Location {L_BIT=1, R_BIT=2, T_BIT=4, B_BIT=8, IN_BIT=16, }; int sphere_location(const T::Sphere_point_3 &sp, - T::Key locate_point, T::Key sphere) const ; - bool behind_arc(const T::Sphere_point_3 &ep, T::Key ind, + bool behind_arc(const T::Sphere_point_3 &ep, Sds::Curve arc, int location) const; void point_sphere_orientation(const T::Sphere_point_3 &time, - T::Key point, T::Key sphere, std::vector &locations) const; bool locate_point_check_face(const T::Sphere_point_3 &z, - Face_const_handle it, - T::Key ind, + Face_const_handle it, std::vector &locations) const ; bool locate_point_check_face_arcs(const T::Sphere_point_3 &z, - T::Key ind, Face_const_handle f, std::vector &locations) const ; bool locate_point_check_face_vertices(const T::Sphere_point_3 &ep, - T::Key index, Face_const_handle it) const; CGAL::Comparison_result debug_rule_shoot_answer(const T::Sphere_point_3 &z, @@ -413,9 +412,7 @@ struct Slice { CGAL::Comparison_result rule_shoot_compare_SS(const T::Sphere_point_3 &z, Sds::Curve srule, - Sds::Curve arc0, - Sds::Point pt, - Sds::Curve arc1) const ; + Sds::Point pt) const ; // return comparison of separator to intersection point on the C coordinate @@ -430,10 +427,10 @@ struct Slice { Sds::Curve rule, Sds::Point pt, CGAL::Comparison_result &ret) const ; - bool rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 &z, + /*bool rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 &z, Sds::Curve rule, Sds::Curve a, - CGAL::Comparison_result &ret) const; + CGAL::Comparison_result &ret) const;*/ @@ -441,7 +438,7 @@ struct Slice { // Debug Functions --------------------------------------------------- - Face_handle locate_point(const T::Sphere_point_3 &ep); + //Face_handle locate_point(const T::Sphere_point_3 &ep); Halfedge_handle shoot_rule(const T::Sphere_point_3 &source, Face_handle f, diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice_data_structure.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice_data_structure.h index bf3e90b502f..5298728c7e9 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice_data_structure.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice_data_structure.h @@ -196,7 +196,9 @@ public: void exchange_vertices(Halfedge_handle h, Halfedge_handle p); - Face_handle intersect(Halfedge_handle ha, Halfedge_handle hb); + std::pair + intersect(Halfedge_handle ha, Halfedge_handle hb, + Halfedge_handle hc, Halfedge_handle hd); std::pair unintersect(Face_handle f); @@ -252,6 +254,11 @@ public: // a halfedge on the rule pointing to the extremal vertex Halfedge_handle rule_halfedge(Curve::Key k, int i) const; + // a halfedge on the circle pointing to the extremal vertex + Halfedge_handle extremum_halfedge(Curve::Key k, int i) const; + + void set_extremum_halfedge(Curve::Key k, int i, Halfedge_handle h); + template void find_halfedges(Curve c, Out o) const { diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_line_intersection.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_line_intersection.h index 0ead582a819..7442ae4a009 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_line_intersection.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_line_intersection.h @@ -177,7 +177,7 @@ public: } CGAL::Comparison_result compare_on_line(const Point_3 &pt) const { - CGAL_precondition(line().has_on(pt)); + CGAL_exactness_precondition(line().has_on(pt)); for (unsigned int i=0; i< 3; ++i){ CGAL::Comparison_result c= compare(pt,Coordinate_index(i)); if (c != CGAL::EQUAL) { diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_traits_3.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_traits_3.h index 0711bff0e6c..be396354441 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_traits_3.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_traits_3.h @@ -139,11 +139,15 @@ struct Arrangement_of_spheres_traits_3 { bool sphere_intersects_rule(Key sphere, Key rule_sphere, Coordinate_index C) const; - bool sphere_intersects_sweep(Key sphere, Sphere_point_3 ep) const; + bool sphere_intersects_rule(Key sphere, const Sphere_point_3& sp, + Coordinate_index C) const; + + bool sphere_intersects_sweep(Key sphere, const Sphere_point_3& ep) const; + CGAL::Comparison_result compare_equipower_point_to_rule(Key sphere0, Key sphere1, - Key rule_sphere, + const Sphere_point_3& sp, Coordinate_index C) const; CGAL::Sign sign_of_separating_plane_normal_c(Key sphere0, Key sphere1, @@ -157,8 +161,8 @@ struct Arrangement_of_spheres_traits_3 { CGAL::Oriented_side oriented_side_of_equipower_plane(Key sphere_0, Key sphere_1, const Sphere_point_3 &s) const; - CGAL::Oriented_side oriented_side_of_center_plane(Key sphere_0, Key sphere_1, - Key sphere_center) const ; + CGAL::Oriented_side oriented_side_of_separating_plane(Key sphere_0, Key sphere_1, + const Sphere_point_3& sp) const ; CGAL::Comparison_result compare_sphere_centers_c(Key a, Key b, Coordinate_index C) const; @@ -170,20 +174,32 @@ struct Arrangement_of_spheres_traits_3 { const Sphere_point_3& d, Coordinate_index C) const; - CGAL::Bounded_side bounded_side_of_sphere_on_equipower_plane_rule_line(Key sphere0, - Key Sphere1, - Key rule, - Coordinate_index C, - const Sphere_point_3 &sp) const; + // return true if the interval defined by the sphere in C contains d.C + bool is_over_circle_c(Key c, const Sphere_point_3& d, + Coordinate_index C) const; + + CGAL::Comparison_result compare_sphere_sphere_at_sweep(Key sphere0, + Key Sphere1, + const Sphere_point_3 &sw, + const Sphere_point_3 &ep, + Coordinate_index C) const; + + CGAL::Bounded_side bounded_side_of_sphere_projected(Key sphere, + const Sphere_point_3 &ep, + Key plane, + Coordinate_index C) const; CGAL::Bounded_side bounded_side_of_sphere(Key sphere, - Key sphere_x, - Key sphere_y, + const Sphere_point_3 &z) const; + CGAL::Bounded_side bounded_side_of_sphere(Key sphere, + Key x, + Key y, const Sphere_point_3 &z) const; CGAL::Comparison_result compare_depths(const Sphere_point_3 &a, const Sphere_point_3 &b) const; - + CGAL::Oriented_side oriented_side(const Plane_3 &p, + const Sphere_point_3 &pt) const; private: // HDS hds_; diff --git a/Arrangement_of_spheres_3/src/CGAL/Arrangement_of_spheres_traits_3.C b/Arrangement_of_spheres_3/src/CGAL/Arrangement_of_spheres_traits_3.C index a2dd3895d50..d36e129f3bc 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Arrangement_of_spheres_traits_3.C +++ b/Arrangement_of_spheres_3/src/CGAL/Arrangement_of_spheres_traits_3.C @@ -23,7 +23,7 @@ bool Arrangement_of_spheres_traits_3::intersects(Key a, Point_3 eqpt= table_->equipower_point(a,b); std::cout << eqpt << std::endl; CGAL::Bounded_side bs=table_->sphere(a).bounded_side(eqpt); - CGAL::Bounded_side obs=table_->sphere(b).bounded_side(eqpt); + //CGAL::Bounded_side obs=table_->sphere(b).bounded_side(eqpt); return bs != CGAL::ON_UNBOUNDED_SIDE; } catch (Sphere_3_table::Equal_centers_exception) { std::cout << "Equal centers for " << a << " and " << b << std::endl; @@ -56,7 +56,7 @@ bool Arrangement_of_spheres_traits_3::intersects(Key a, bool Arrangement_of_spheres_traits_3::sphere_intersects_sweep(Key sphere, - Sphere_point_3 ep) const { + const Sphere_point_3& ep) const { Event_pair be= sphere_events(sphere); if (be.first.compare(ep, sweep_coordinate()) != CGAL::LARGER && be.second.compare(ep, sweep_coordinate()) != CGAL::SMALLER) return true; @@ -77,6 +77,24 @@ bool Arrangement_of_spheres_traits_3::sphere_intersects_rule(Key s, } +bool Arrangement_of_spheres_traits_3::sphere_intersects_rule(Key s, + const Sphere_point_3 &sp, + Coordinate_index C) const { + FT v[3]={0,0,0}; + CGAL::Comparison_result cr=sp.compare(table_->center(s), C); + if (cr == CGAL::SMALLER) { + v[C.index()]=-1; + } else if (cr == CGAL::LARGER) { + v[C.index()]=1; + } else { + return true; + } + Sphere_point_3 ep(table_->sphere(s), Line_3(table_->center(s), + Vector_3(v[0], v[1], v[2]))); + return ep.compare(sp, C) != cr; +} + + CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_sphere_center_c(Key a, const Sphere_point_3& d, @@ -95,12 +113,11 @@ Arrangement_of_spheres_traits_3::compare_sphere_center_c(Key a, CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_equipower_point_to_rule(Key a, Key b, - Key c, + const Sphere_point_3 &sp, Coordinate_index C) const{ // NOTE redo this with computing it directly Point_3 pt= table_->equipower_point(a,b); - return CGAL::compare(pt[C.index()], table_->center(c)[C.index()]); - + return CGAL::Comparison_result(-sp.compare(pt, C)); } @@ -129,19 +146,16 @@ CGAL::Sign Arrangement_of_spheres_traits_3::sign_of_equipower_plane_normal_c(Key return sn; } - - -CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_plane(Key a, Key b, - const Sphere_point_3 &s) const { - CGAL_assertion(s.is_valid()); - Plane_3 p= table_->equipower_plane(a,b); +CGAL::Oriented_side +Arrangement_of_spheres_traits_3::oriented_side(const Plane_3 &p, + const Sphere_point_3 &s) const{ CGAL::Object o= di_(p, s.line()); Point_3 pt; Line_3 l; if (CGAL::assign(pt, o)){ CGAL::Comparison_result c= s.compare_on_line(pt); - FT d= p.orthogonal_vector() * s.line().to_vector(); - if (d > 0) return CGAL::enum_cast(c); + CGAL::Sign sn= CGAL::sign(p.orthogonal_vector() * s.line().to_vector()); + if (sn ==CGAL::POSITIVE) return CGAL::enum_cast(c); else return CGAL::enum_cast(-c); } else if (CGAL::assign( l, o)){ return CGAL::ON_ORIENTED_BOUNDARY; @@ -150,15 +164,19 @@ CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_ } } +CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_plane(Key a, Key b, + const Sphere_point_3 &s) const { + CGAL_assertion(s.is_valid()); + Plane_3 p= table_->equipower_plane(a,b); + return oriented_side(p, s); +} -CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_center_plane(Key a, Key b, - Key sphere_center) const { - Vector_3 d(table_->center(b)-table_->center(a)); - Line_2 l(Point_2(table_->center(a)[plane_coordinate(0).index()], - table_->center(a)[plane_coordinate(1).index()]), - Vector_2(d[plane_coordinate(0).index()], d[plane_coordinate(1).index()])); - return l.oriented_side(Point_2(table_->center(sphere_center)[plane_coordinate(0).index()], - table_->center(sphere_center)[plane_coordinate(1).index()])); + +CGAL::Oriented_side +Arrangement_of_spheres_traits_3::oriented_side_of_separating_plane(Key a, Key b, + const Sphere_point_3 &ep) const { + Plane_3 sp= table_->separating_plane(a,b); + return oriented_side(sp, ep); } CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_sphere_centers_c(Key a, Key b, Coordinate_index C) const { @@ -167,56 +185,192 @@ CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_sphere_centers_ +CGAL::Comparison_result +Arrangement_of_spheres_traits_3::compare_sphere_sphere_at_sweep(Key sphere0, + Key sphere1, + const Sphere_point_3 &sw, + const Sphere_point_3 &ep, + Coordinate_index C) const { + { + std::cout << "Compare SS at sweep " << sphere0 << " " << sphere1 + << " sweep " << sw << " line " << ep.line() << " coord " << C + << std::endl; + + bool i0= sphere_intersects_rule(sphere0, ep, C); + bool i1= sphere_intersects_rule(sphere1, ep, C); + + + //std::cout << "Comparisons are " << hi << " " << hp << " and " << ohi << " " << ohp << std::endl; + if (!i0) { + CGAL::Comparison_result sc0= compare_sphere_center_c(sphere0, ep, C); + std::cout << "Returning " << sc0 << " since sphere 0 misses" << std::endl; + return sc0; + } else if (!i1) { + CGAL::Comparison_result sc1= compare_sphere_center_c(sphere1, ep, C); + std::cout << "Returning " << sc1 << " since sphere 0 misses" << std::endl; + return sc1; + } + } + + + Plane_3 eqp= table_->equipower_plane(sphere0, sphere1); + std::cout << "EQP is " << eqp << std::endl; + + FT v[3]={0,0,0}; + v[other_plane_coordinate(C).index()]=1; + + Plane_3 rule_plane(ep.line().point(), + ep.line().point()+ Vector_3(v[0], v[1], v[2]), + ep.line().point()+ ep.line().to_vector()); + if (rule_plane.orthogonal_vector()[C.index()] < 0) { + rule_plane=rule_plane.opposite(); + } + std::cout << "rule_plane is " << rule_plane << std::endl; + CGAL_assertion(rule_plane != Plane_3()); + Line_3 l; + CGAL::Object o= di_(eqp, rule_plane); + if (!CGAL::assign(l,o)) { + CGAL::Comparison_result cr= compare_equipower_point_to_rule(sphere0, sphere1, + ep, C); + std::cout << "Returning " << cr << " since sphere the planes miss" << std::endl; + return cr; + } + + std::cout << "The intersection line is " << l << std::endl; + + + + Sphere_point_3 switch_point_0(table_->sphere(sphere0), l); + + if (!switch_point_0.is_valid()){ + + CGAL::Comparison_result cr= compare_equipower_point_to_rule(sphere0,sphere1, + ep, C); + std::cout << "Returning " << cr << " line misses" << std::endl; + return cr; + } + + //if (switch_point_0 > switch_point_1) { + // std::swap(switch_point_0, switch_point_1); + //} + Plane_3 sepp= table_->separating_plane(sphere0, sphere1); + std::cout << "separating_plane is " << sepp << std::endl; + + Sphere_point_3 frontp= intersection_2_events(sphere0, sphere1).first; + std::cout << "The front point is " << frontp << std::endl; + CGAL::Oriented_side fpos= oriented_side(rule_plane, frontp); + bool flip=false; + if (fpos == CGAL::ON_ORIENTED_BOUNDARY) { + std::cout << "Front point is on boundary using normals "; + fpos = CGAL::Oriented_side(CGAL::sign(sepp.orthogonal_vector() + * rule_plane.orthogonal_vector())); + std::cout << "got " << fpos << std::endl; + } + + std::cout << "Oriented_side is " << fpos << std::endl; + CGAL::Comparison_result spcr0= switch_point_0.compare(sw, sweep_coordinate()); + + Sphere_point_3 switch_point_1(table_->sphere(sphere0), l.opposite()); + std::cout << "switch points are " << switch_point_0 + << " and " << switch_point_1 << std::endl; + CGAL::Comparison_result spcr1= switch_point_1.compare(sw, sweep_coordinate()); + std::cout << "spcr is " << spcr0 << spcr1 << std::endl; + + + CGAL::Oriented_side os_0= oriented_side(sepp, switch_point_0); + CGAL::Oriented_side os_1= oriented_side(sepp, switch_point_1); + std::cout << "os are " << os_0 << os_1 << std::endl; + + if (os_0 == CGAL::ON_POSITIVE_SIDE && os_1 == CGAL::ON_POSITIVE_SIDE) { + flip = (spcr0 != spcr1); + } else if (os_0 == CGAL::ON_POSITIVE_SIDE) { + flip = (spcr0 == CGAL::SMALLER); + } else if (os_1 == CGAL::ON_POSITIVE_SIDE) { + flip = (spcr1 == CGAL::LARGER); + } + + if (!flip) { + return CGAL::Comparison_result(fpos); + } else { + return CGAL::Comparison_result(-fpos); + } + } + + CGAL::Bounded_side -Arrangement_of_spheres_traits_3::bounded_side_of_sphere_on_equipower_plane_rule_line(Key sphere0, - Key sphere1, - Key rule, - Coordinate_index C, - const Sphere_point_3& sp) const { - - Plane_3 eqp= table_->equipower_plane(sphere0, sphere1); - FT v[3], p[3]; - v[C.index()]=1; - v[other_plane_coordinate(C).index()]=0; - v[sweep_coordinate().index()] = 0; - p[plane_coordinate(0).index()]=table_->center(rule)[plane_coordinate(0).index()]; - p[plane_coordinate(1).index()]=table_->center(rule)[plane_coordinate(1).index()]; - p[sweep_coordinate().index()]= 0; - Plane_3 pp(Point_3(p[0],p[1],p[2]), - Vector_3(v[0], v[1], v[2])); - Line_3 l; - CGAL::Object o= di_(eqp, pp); - CGAL_assertion(CGAL::assign(l,o)); - CGAL::assign(l,o); - std::cout << "The intersection line is " << l << std::endl; - //CGAL::Comparison_result dir= CGAL::compare(ep.simple_coordinate(1-C), spheres_[h->curve().key()].table_->center(1-C));*/ - Sphere_point_3 p0(table_->sphere(sphere0), l); +Arrangement_of_spheres_traits_3 +::bounded_side_of_sphere_projected(Key sphere, + const Sphere_point_3 &ep, + Key plane, + Coordinate_index C) const{ + FT p[3],v[3]; + for (unsigned int i=0; i< 3; ++i) { + p[i]= ep.line().point()[i]; + v[i]= ep.line().to_vector()[i]; + } + v[other_plane_coordinate(C).index()] = 0; + p[other_plane_coordinate(C).index()] = table_->center(plane)[other_plane_coordinate(C).index()]; - - if (!p0.is_valid()){ + Line_3 l(Point_3(p[0],p[1], p[2]),Vector_3(v[0], v[1], v[2])); + std::cout << "Line is " << l << std::endl; + Sphere_point_3 p0(table_->sphere(sphere), l); + if (!p0.is_valid()) { + std::cout << "Line misses" << std::endl; return CGAL::ON_UNBOUNDED_SIDE; } else { - Sphere_point_3 p1(table_->sphere(sphere0), l.opposite()); - - if (compare_depths(p0, p1)== CGAL::LARGER) { - std::swap(p0, p1); - } - //CGAL_assertion(p0 < p1); - if (compare_depths(sp,p0)==CGAL::EQUAL - || compare_depths(sp, p1) == CGAL::EQUAL) { + Sphere_point_3 p1(table_->sphere(sphere), l.opposite()); + std::cout << "Points are " << p0 << " " << p1 << std::endl; + CGAL::Comparison_result c0= compare_depths(ep, p0), c1= compare_depths(ep, p1); + std::cout << "Depths are " << c0 << " " << c1 << std::endl; + if (c0==CGAL::EQUAL || c1 == CGAL::EQUAL) { return CGAL::ON_BOUNDARY; + } else if (c0 != c1) { + return CGAL::ON_BOUNDED_SIDE; + } else { + return CGAL::ON_UNBOUNDED_SIDE; } - if (compare_depths(sp, p0) == CGAL::LARGER - && compare_depths(sp, p1) == CGAL::SMALLER) return CGAL::ON_BOUNDED_SIDE; - else return CGAL::ON_UNBOUNDED_SIDE; } } - // compare along C where the intersection point of the two rules // (defined by Coord C of ep) is relative to the sphere of the arc +CGAL::Bounded_side +Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s, + const Sphere_point_3 &z) const{ + /*FT p[3]; + p[plane_coordinate(0).index()]= table_->center(x)[plane_coordinate(0).index()]; + p[plane_coordinate(1).index()]= table_->center(y)[plane_coordinate(1).index()]; + p[sweep_coordinate().index()]=0; + FT v[3]; + v[plane_coordinate(0).index()]=0; + v[plane_coordinate(1).index()]=0; + v[sweep_coordinate().index()]=1;*/ + Line_3 l= z.line();//(Point_3(p[0], p[1], p[2]), Vector_3(v[0],v[1],v[2])); + + Sphere_point_3 p0(table_->sphere(s), l); + + if (!p0.is_valid()) { + return CGAL::ON_UNBOUNDED_SIDE; + } else { + Sphere_point_3 p1(table_->sphere(s), + l.opposite()); + CGAL::Comparison_result c0= compare_depths(p0, z); + CGAL::Comparison_result c1= compare_depths(p1, z); + if (c0 == CGAL::EQUAL || c1== CGAL::EQUAL) { + //std::cout << "Equal" << std::endl; + return CGAL::ON_BOUNDARY; + } + //CGAL_assertion(p0 <= p1); + if (c0== CGAL::SMALLER && c1 == CGAL::LARGER){ + return CGAL::ON_BOUNDED_SIDE; + } else { + return CGAL::ON_UNBOUNDED_SIDE; + } + } +} + CGAL::Bounded_side Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s, Key x, @@ -230,8 +384,6 @@ Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s, v[plane_coordinate(0).index()]=0; v[plane_coordinate(1).index()]=0; v[sweep_coordinate().index()]=1; - /*std::cout << "Line through " << CGAL::to_double(p[0]) - << " " << CGAL::to_double(p[1]) << std::endl;*/ Line_3 l(Point_3(p[0], p[1], p[2]), Vector_3(v[0],v[1],v[2])); Sphere_point_3 p0(table_->sphere(s), l); @@ -258,6 +410,32 @@ Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s, +bool +Arrangement_of_spheres_traits_3::is_over_circle_c(Key s, + const Sphere_point_3 &z, + Coordinate_index C) const { + Line_3 l= z.line(); + FT p[3], v[3]; + for (unsigned int i=0; i< 3; ++i) { + p[i]= l.point()[i]; + v[i]= l.to_vector()[i]; + } + Coordinate_index oc= other_plane_coordinate(C); + p[oc.index()]= table_->center(s)[oc.index()]; + v[oc.index()]= 0; + Line_3 nl=Line_3(Point_3(p[0], p[1], p[2]),Vector_3(v[0], v[1], v[2])); + Sphere_point_3 sp(table_->sphere(s), nl); + if (!sp.is_valid()) return false; + else { + Sphere_point_3 spo(table_->sphere(s), nl.opposite()); + CGAL::Comparison_result c= sp.compare(z, sweep_coordinate()); + CGAL::Comparison_result co= spo.compare(z, sweep_coordinate()); + if (co == c) return true; + else return false; + } +} + + /* diff --git a/Arrangement_of_spheres_3/src/CGAL/Combinatorial_stuff.C b/Arrangement_of_spheres_3/src/CGAL/Combinatorial_stuff.C index e621cbc16ec..722e68d696c 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Combinatorial_stuff.C +++ b/Arrangement_of_spheres_3/src/CGAL/Combinatorial_stuff.C @@ -344,6 +344,7 @@ Combinatorial_curve::constant_coordinate() const { r= plane_coordinate(1); else r= plane_coordinate(0); //if (!is_finite()) r=1-r; + //std::cout << "Constant coordinate of " << *this << " is " << r << std::endl; return r; } diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_arrangement.C b/Arrangement_of_spheres_3/src/CGAL/Slice_arrangement.C index ac5d449e4b3..45b0427b6f7 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_arrangement.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_arrangement.C @@ -544,11 +544,13 @@ Slice_arrangement::Point Slice_arrangement::point(CArr::Vertex_const_handle h) c } Slice_arrangement::Curve Slice_arrangement::curve(CArr::Halfedge_const_handle h) const { - CGAL_assertion(std::distance(arr_.originating_curves_begin(h), - arr_.originating_curves_end(h))==1); - CGAL_assertion(map_.find(arr_.originating_curves_begin(h)) != map_.end()); - return map_.find(arr_.originating_curves_begin(h))->second; + Slice_arrangement::Curve c= map_.find(arr_.originating_curves_begin(h))->second; + if (std::distance(arr_.originating_curves_begin(h), + arr_.originating_curves_end(h))!=1) { + CGAL_assertion(c.is_rule()); + } + return c; //return boost::apply_visitor(lookup_curve_data_, h->curve()); /*if (Circular_k::Circular_arc_2& a=boost::get(h->curve())) { diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_data_structure.C b/Arrangement_of_spheres_3/src/CGAL/Slice_data_structure.C index 92cfb673b08..ac02410b8cc 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_data_structure.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_data_structure.C @@ -36,6 +36,24 @@ Slice_data_structure::rule_halfedge(Curve::Key k, int i) const { return h; } +Slice_data_structure::Halfedge_handle +Slice_data_structure::extremum_halfedge(Curve::Key k, int i) const { + Halfedge_handle h= halfedges_[k.input_index()][i]; + CGAL_assertion(h->vertex()->point().is_sphere_extremum()); + CGAL_assertion(h->curve().key() == k); + CGAL_assertion(!h->curve().is_rule()); + return h; +} + +void +Slice_data_structure::set_extremum_halfedge(Curve::Key k, int i, + Halfedge_handle h){ + CGAL_precondition(h->curve().key() ==k); + CGAL_precondition(h->curve().is_arc()); + CGAL_precondition(h->vertex()->point().is_sphere_extremum()); + halfedges_[k.input_index()][i]; +} + Slice_data_structure::Halfedge_handle Slice_data_structure::next_edge_on_curve(Halfedge_handle h) const { if (h->curve().is_rule()) { @@ -271,7 +289,7 @@ Slice_data_structure::exchange_spheres(Curve::Key k, Curve::Key l) { if (kh->vertex()->point().is_sphere_extremum()) { // relabel rule Halfedge_handle r= kh->opposite()->prev()->opposite(); - kes[r->curve().rule_index()]= r; + //kes[r->curve().rule_index()]= r; CGAL_assertion(r->curve().is_rule()); if (r->curve().key() == k) { Curve nr= r->curve(); @@ -646,6 +664,7 @@ bool Slice_data_structure::has_vertex(Face_const_handle fh, Halfedge_const_handle h= vh->halfedge(); do { if (h->face() == fh) return true; + //if (h->opposite()->face() == fh) return true; h= h->opposite()->prev(); CGAL_assertion(h->vertex() == vh); } while (h != vh->halfedge()); @@ -1203,8 +1222,10 @@ void Slice_data_structure::exchange_vertices(Halfedge_handle h, -Slice_data_structure::Face_handle -Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb) { +std::pair +Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb, + Halfedge_handle hc, Halfedge_handle hd) { CGAL_precondition(ha->face() == hb->face()); Halfedge_handle ha_prev= ha->prev(); Halfedge_handle ha_op_next= ha->opposite()->next(); @@ -1302,7 +1323,8 @@ Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb) { write_face(hb->face()->halfedge(), std::cout); std::cout << std::endl; - return fn; + CGAL_assertion(0); + return std::pair();; } diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_kds.C b/Arrangement_of_spheres_3/src/CGAL/Slice_kds.C index 999d7fbf71f..4113fbce993 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_kds.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_kds.C @@ -212,43 +212,55 @@ void Slice::process_one_sphere_event(T::Key k) { } } -void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { - std::cout << "Two sphere kds " << k << " " << l << " " << f << std::endl; +void Slice::process_equal_centers_two_sphere_event(T::Key k, T::Key l, bool f) { Intersection_2 i2(k,l); - if (!f) { - intersections_2_.erase(i2); - } else { - intersections_2_[i2].first=Event_key(); + if (t_.compare_sphere_centers_c(k,l, sweep_coordinate()) == CGAL::LARGER) { + std::swap(k,l); } - // check if centers align - CGAL::Comparison_result c[2]; - c[0]= t_.compare_sphere_centers_c(k,l, plane_coordinate(0)); - c[1] t_.compare_sphere_centers_c(k,l, plane_coordinate(1)); - if (c[0]== CGAL::EQUAL && c[1] == CGAL::EQUAL) { - if (t_.compare_sphere_centers_c(k,l, sweep_coordinate()) == CGAL::LARGER) { - std::swap(k,l); - } - sds_.exchange_spheres(k,l); - sim_->delete_event(intersections_2_[i2].second); - intersections_2_.erase(i2); - - for (unsigned int i=0; i< 4; ++i) { - Halfedge_handle rh= sds_.rule_halfedge(k, i); - // now outside - clean_edge(rh); - check_edge_collapse(rh); - } - for (unsigned int i=0; i< 4; ++i) { - Halfedge_handle rh= sds_.rule_halfedge(l, i); - // now outside - clean_edge(rh); - if (!(rh->vertex()->point().is_sphere_rule() - && rh->vertex()->point().sphere_key() == k)) { + sds_.exchange_spheres(k,l); + sim_->delete_event(intersections_2_[i2].second); + intersections_2_.erase(i2); + { + Halfedge_handle kc= sds_.a_halfedge(k), ke=kc; + do { + kc->curve().set_key(l); + kc= sds_.next_edge_on_curve(kc); + } while (kc != ke); + } + { + Halfedge_handle lc= sds_.a_halfedge(l), le=lc; + do { + lc->curve().set_key(k); + lc= sds_.next_edge_on_curve(lc); + } while (lc != le); + } + Halfedge_handle keh[4]; + for (unsigned int i=0; i< 4; ++i) { + Halfedge_handle rh= sds_.rule_halfedge(k, i); + keh[i]= rh; + sds_.set_extremum_halfedge(k,i,sds_.extremum_halfedge(l,i)); + // now outside + clean_edge(rh); + check_edge_collapse(rh); + } + for (unsigned int i=0; i< 4; ++i) { + Halfedge_handle rh= sds_.rule_halfedge(l, i); + // now outside + sds_.set_extremum_halfedge(l,i,keh[i]); + clean_edge(rh); + if (!(rh->vertex()->point().is_sphere_rule() + && rh->vertex()->point().sphere_key() == k)) { check_edge_collapse(rh); - } } - } else if (c[0]== CGAL::EQUAL || c[1]== CGAL::EQUAL) { - // figure out which extremums + } +} + + + +void Slice::process_aligned_centers_two_sphere_event(T::Key k, T::Key l, bool f, + CGAL::Comparison_result c[2]) { + // figure out which extremums + Intersection_2 i2(k,l); const int equal_dir= (c[0]== CGAL::EQUAL)? 0:1; const int other_dir= 1-equal_dir; // make sure that k is smaller @@ -272,25 +284,38 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { // name them for who connects later Halfedge_handle rule_vertex_l; - if (extr_k->next()->key() != k) { + if (extr_k->next()->curve().key() != k) { clean_edge(extr_k->next()); std::pair hr - = sds_.remove_rule(extra_k->next()); + = sds_.remove_rule(extr_k->next()); rule_vertex_l= hr.first; CGAL_assertion(hr.second == extr_k); } Halfedge_handle rule_vertex_k; - if (extr_l->next()->key() != l) { + if (extr_l->next()->curve().key() != l) { clean_edge(extr_l->next()); std::pair hr - = sds_.remove_rule(extra_l->next()); + = sds_.remove_rule(extr_l->next()); rule_vertex_k= hr.first; CGAL_assertion(hr.second == extr_l); } // maybe replace by two pinch operations - std::pair hp= - intersect_2(extr_l->next(), extr_l->opposite(), - extr_k->next(), extr_k->opposite()); + clean_edge(extr_l->opposite()); + clean_edge(extr_k->opposite()); + CGAL_assertion(0); + std::pair hp; + /*= + sds_.intersect_2(extr_l->next(), extr_l->opposite(), + extr_k->next(), extr_k->opposite());*/ + // clean edges and create certificates + CGAL_postcondition(hp.first->face() == hp.second->face()); + CGAL_postcondition(hp.first->prev()->prev()== hp.second); + CGAL_postcondition(hp.second->prev()->prev()== hp.first); + + check_edge_collapse(hp.first->opposite()->prev()); + check_edge_collapse(hp.first->opposite()->prev()->opposite()->prev()); + check_edge_collapse(hp.second->opposite()->prev()); + check_edge_collapse(hp.second->opposite()->prev()->opposite()->prev()); CGAL_assertion(hp.first->curve().key() == k); CGAL_assertion(hp.second->curve().key() == l); @@ -310,11 +335,11 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { Sds::Curve::make_rule(l, rl)); // set halfedges for spheres - sds_.set_extremum_edge(k, hp.first->prev()); - sds_.set_extremum_edge(l, hp.second->prev()); + sds_.set_extremum_halfedge(k, rk, hp.first->prev()); + sds_.set_extremum_halfedge(l, rl, hp.second->prev()); } else { // first figure out which case we are in - const T::Event_point_3 ep& = sim_->current_time(); + const T::Event_point_3 &ep = sim_->current_time(); CGAL::Comparison_result c[2]; c[0]= t_.compare_sphere_center_c(k, ep, T::Coordinate_index(other_dir)); c[1]= t_.compare_sphere_center_c(l, ep, T::Coordinate_index(other_dir)); @@ -327,6 +352,26 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { // then handle hard case // handle easy case by accident } +} + + +void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { +#if 0 + std::cout << "Two sphere kds " << k << " " << l << " " << f << std::endl; + Intersection_2 i2(k,l); + if (!f) { + intersections_2_.erase(i2); + } else { + intersections_2_[i2].first=Event_key(); + } + // check if centers align + CGAL::Comparison_result c[2]; + c[0]= t_.compare_sphere_centers_c(k,l, plane_coordinate(0)); + c[1] t_.compare_sphere_centers_c(k,l, plane_coordinate(1)); + if (c[0]== CGAL::EQUAL && c[1] == CGAL::EQUAL) { + process_equal_centers_two_sphere_event(k,l,f); + } else if (c[0]== CGAL::EQUAL || c[1]== CGAL::EQUAL) { + process_aligned_centers_two_sphere_event(k,l,f, c); } else if (f) { typedef std::pair EP; @@ -380,6 +425,7 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) { //CGAL_assertion(!shared_faces.empty()); unintersect_spheres(sim_->current_time(), f); } +#endif } void Slice::process_three_sphere_event(T::Key k, T::Key l, T::Key m, bool f) { diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_modifiers.C b/Arrangement_of_spheres_3/src/CGAL/Slice_modifiers.C index 982dff14a22..3d2ac2448b4 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_modifiers.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_modifiers.C @@ -389,7 +389,7 @@ Slice::Face_handle Slice::insert_sphere(const T::Sphere_point_3 &ep, // the halfedge points to the vertex and is on this face Vertex_handle vhs[4]; try { - f= locate_point(ep, k); + f= locate_point(ep); } catch (On_edge_exception e) { std::cout << "Point hit edge "; sds_.write( e.halfedge_handle(), std::cout) << std::endl; diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_point_location.C b/Arrangement_of_spheres_3/src/CGAL/Slice_point_location.C index 36ddd5a83ee..e655f7b6c84 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_point_location.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_point_location.C @@ -15,7 +15,6 @@ bool Slice::locate_point_check_face(const T::Sphere_point_3 &z, Face_const_handle it, - T::Key index, std::vector &locations) const { Halfedge_const_handle h= it->halfedge(); //bool finite=false; @@ -26,12 +25,12 @@ bool Slice::locate_point_check_face(const T::Sphere_point_3 &z, if (h->curve().is_right() && h->curve().is_inside()) return false; } else { T::Key sphere= h->curve().key(); - point_sphere_orientation(z, index, sphere, locations); + point_sphere_orientation(z, sphere, locations); //DPRINT(std::cout << "Testing " << h->curve() <curve().is_compatible_location(locations[sphere.input_index()])) { //DPRINT(std::cout << "Nixed by edge " << h->curve() << std::endl); return false; - } + } } //finite=true; //} else { @@ -57,7 +56,6 @@ bool Slice::locate_point_check_face(const T::Sphere_point_3 &z, // cache checked arcs bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep, - T::Key ind, Face_const_handle f, std::vector &locations) const { std::set check_arcs; @@ -67,7 +65,7 @@ bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep, && check_arcs.find(h->curve().key()) == check_arcs.end()){ std::cout << "Arc test " << ep << " on " << h->curve() << std::endl; check_arcs.insert(h->curve().key()); - bool ba=behind_arc(ep, ind, h->curve(), + bool ba=behind_arc(ep, h->curve(), locations[h->curve().key().input_index()]); if (ba) { std::cout << "Point is behind arc " << std::endl; @@ -84,8 +82,8 @@ bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep, bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep, - T::Key index, Face_const_handle it) const { + Halfedge_const_handle h= it->halfedge(); do { if (h->vertex()->point().is_sphere_sphere() @@ -94,17 +92,17 @@ bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep, // NOTE what about degeneracies? not sure if I need to handle // them here Sds::Point npt= h->vertex()->point(); - if (t_.oriented_side_of_center_plane(npt.sphere_key(0), - npt.sphere_key(1), - index) == CGAL::ON_NEGATIVE_SIDE) { - CGAL_assertion(0); + if (t_.oriented_side_of_separating_plane(npt.sphere_key(0), + npt.sphere_key(1), + ep) == CGAL::ON_NEGATIVE_SIDE) { + //CGAL_assertion(0); std::cout << "Face nixed by vertex " << npt << std::endl; return false; } } h= h->next(); } while (h != it->halfedge()); - + return true; } @@ -126,18 +124,17 @@ bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep, -Slice::Face_handle Slice::locate_point(const T::Sphere_point_3 & ep) { +/*Slice::Face_handle Slice::locate_point(const T::Sphere_point_3 & ep) { t_.set_temp_sphere(ep.sphere()); return locate_point(ep, T::Key::temp_key()); -} + }*/ Slice::Face_handle - Slice::locate_point(const T::Sphere_point_3 & ep, - T::Key index) { - return locate_point(sds_.faces_begin(), sds_.faces_end(), ep, index); + Slice::locate_point(const T::Sphere_point_3 & ep) { + return locate_point(sds_.faces_begin(), sds_.faces_end(), ep); } diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_predicates.C b/Arrangement_of_spheres_3/src/CGAL/Slice_predicates.C index 1d351caa1b1..a52c12bafd6 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_predicates.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_predicates.C @@ -13,12 +13,9 @@ int Slice::sphere_location(const T::Sphere_point_3& sp, - T::Key locate_point, T::Key s) const { int r=0; CGAL::Bounded_side bs=t_.bounded_side_of_sphere(s, - locate_point, - locate_point, sp); /*T::Event_point_3 ep=sphere_start(locate_point); @@ -37,17 +34,15 @@ int Slice::sphere_location(const T::Sphere_point_3& sp, } else { r |= Sds::Curve::lOUT_BIT; } - - /*T::Plane_3 lrp(center(s), T::Vector_3(1, 0, 0)); - T::Plane_3 tbp(center(s), T::Vector_3(0, 1, 0)); - CGAL::Comparison_result xo= oriented_side(lrp, ep); - CGAL::Oriented_side yo= oriented_side(tbp, ep);*/ - //CGAL::Comparison_result xo= compare_event_point_to_rule(ep, s, 0); - //CGAL::Comparison_result yo= compare_event_point_to_rule(ep, s, 1); - CGAL::Comparison_result xo= t_.compare_sphere_centers_c(s, locate_point, + + /* + NOTE can make it more specific for intersection point location + */ + + CGAL::Comparison_result xo= t_.compare_sphere_center_c(s, sp, plane_coordinate(0)); - CGAL::Comparison_result yo= t_.compare_sphere_centers_c(s, locate_point, - plane_coordinate(1)); + CGAL::Comparison_result yo= t_.compare_sphere_center_c(s, sp, + plane_coordinate(1)); if (xo != CGAL::LARGER) { r |= Sds::Curve::lR_BIT; } @@ -72,39 +67,19 @@ int Slice::sphere_location(const T::Sphere_point_3& sp, bool Slice::behind_arc(const T::Sphere_point_3 &ep, - T::Key ind, Sds::Curve arc, + Sds::Curve arc, int location) const{ - CGAL_assertion(arc.is_compatible_location(location)); + CGAL_precondition(!arc.is_inside()); + CGAL_precondition(arc.is_compatible_location(location)); T::Coordinate_index C=arc.is_weakly_incompatible(location); - if (C.is_valid()) { - CGAL::Bounded_side bs; - if (C== plane_coordinate(0)) { - bs= t_.bounded_side_of_sphere(arc.key(), - arc.key(), - ind, - ep); - } else { - bs= t_.bounded_side_of_sphere(arc.key(), - ind, - arc.key(), - ep); - } - /*NT v[2]; - v[C]= t_.center_c(arc.key(), C); - v[1-C]= t_.center_c(ind, 1-C); - T::Line_3 l(T::Point_3(v[0], v[1], 0), - T::Vector_3(0,0,1)); - T::Event_point_3 fp(sphere(arc.key()), - l); - if (!fp.is_valid()) return false; - T::Event_point_3 bp(sphere(arc.key()), - l.opposite()); - if (fp <= ep && bp >=ep) { - return true; - } else {return false;} - } else return false;*/ - return bs== CGAL::ON_BOUNDED_SIDE; + if (C.is_valid()) { // we know it is opposite the center + /*CGAL::Comparison_result cr = t_.compare_sphere_center_c(arc.key(), + ep, C); + if (cr == CGAL::SMALLER && arc.is_negative() + || cr == CGAL::LARGER && !arc.is_negative()) return false;*/ + bool b= t_.is_over_circle_c(arc.key(), ep, C); + return b; } else { return false; } @@ -112,13 +87,12 @@ bool Slice::behind_arc(const T::Sphere_point_3 &ep, void Slice::point_sphere_orientation(const T::Sphere_point_3 &time, - T::Key front_point, T::Key sphere, std::vector &locations /*, std::vector &edges*/) const { if (locations[sphere.input_index()]==0){ - locations[sphere.input_index()]=sphere_location(time, front_point, sphere); + locations[sphere.input_index()]=sphere_location(time, sphere); //DPRINT(std::cout << "For sphere " << sphere << " got " << SL::decode(locations[sphere]) << std::endl); } } diff --git a/Arrangement_of_spheres_3/src/CGAL/Slice_shoot_rule.C b/Arrangement_of_spheres_3/src/CGAL/Slice_shoot_rule.C index e274e401ac9..f9fc0c7b0da 100644 --- a/Arrangement_of_spheres_3/src/CGAL/Slice_shoot_rule.C +++ b/Arrangement_of_spheres_3/src/CGAL/Slice_shoot_rule.C @@ -11,93 +11,55 @@ Slice::Halfedge_handle Slice::shoot_rule(const T::Sphere_point_3& source, t_.set_temp_sphere(source.sphere()); return shoot_rule(source, f, Sds::Curve(T::Key::target_key(), Sds::Curve::Part(type))); - } +} + + + Slice::Halfedge_handle Slice::shoot_rule(const T::Sphere_point_3 & source, - Face_handle f, - Sds::Curve rule) { + Face_handle f, + Sds::Curve rule) { bool backwards= (rule.is_top() || rule.is_left()); DPRINT(std::cout << "Rule is " << rule << std::endl); - - Halfedge_handle h= f->halfedge(); - std::vector edges; - - bool has_last=false; - CGAL::Comparison_result last; - if (rule.can_intersect(h->curve()) - && rule.can_intersect(h->prev()->curve())) { - has_last=true; - last= rule_shoot_edge_vertex(source, rule, h->prev()->curve(), - h->prev()->vertex()->point(), - h->curve()); - } else { - DPRINT(std::cout << "Eliminating " << h->prev()->vertex()->point() << std::endl); + + Halfedge_handle end= f->halfedge(); + + while (!rule.can_intersect(end->curve())) { + end=end->next(); } + while (rule.can_intersect(end->curve())) { + end=end->next(); + } + Halfedge_handle h=end; + while (!rule.can_intersect(h->curve())) { + h=h->next(); + } + do { // later skip second one if first is wrong - CGAL::Comparison_result pbs, nbs; - DPRINT(std::cout << "\nTesting edge: " << h->prev()->vertex()->point() << "--" - << h->curve() << "--" - << h->vertex()->point() << std::endl;) - if (rule.can_intersect(h->curve())) { - if (has_last) { - DPRINT(std::cout << "Using last of " << last << std::endl); - pbs=last; - } else if (backwards) { - pbs= CGAL::LARGER; - } else { - pbs=CGAL::SMALLER; - } - /* pbs= shoot_rule_edge_vertex(ep, rule, h->prev()->curve(), - h->prev()->vertex()->point(), - h->curve()); - }*/ - if (rule.can_intersect(h->next()->curve())) { - nbs= rule_shoot_edge_vertex(source, rule, h->curve(), - h->vertex()->point(), - h->next()->curve()); - DPRINT(std::cout << "Results are " << nbs << " " << pbs << std::endl); - has_last=true; - last = nbs; - } else { - has_last=false; - if (backwards) nbs = CGAL::SMALLER; - else nbs = CGAL::LARGER; - } - - if (!backwards && pbs != CGAL::LARGER && nbs != CGAL::SMALLER - || backwards && pbs != CGAL::SMALLER && nbs != CGAL::LARGER) { - //std::cout << "it passes" << std::endl; - edges.push_back(h); - } + DPRINT(sds_.write(h, std::cout) << " is being tested" << std::endl); + if (h->next() == end) { + DPRINT(sds_.write(h, std::cout) << " is returned by default" << std::endl); + return h; } else { - //std::cout << "Filtered" << std::endl; - has_last=false; + CGAL::Comparison_result cr= rule_shoot_edge_vertex(source, rule, h->curve(), + h->vertex()->point(), + h->next()->curve()); + DPRINT(std::cout << "Result is " << cr << std::endl); + if (cr == CGAL::EQUAL) { + DPRINT(std::cout << "On vertex " << h->vertex()->point() << std::endl); + throw On_vertex_exception(h->vertex()); + } + if (backwards && cr == CGAL::SMALLER + || !backwards && cr == CGAL::LARGER) { + DPRINT(sds_.write(h, std::cout) << " is returned " << std::endl); + return h; + } } - - h= h->next(); - CGAL_assertion(h->face()== f); - } while (h != f->halfedge()); + } while (h != end); - - if (edges.size() == 0) { - std::cerr << "No edge found " << std::endl; - return Halfedge_handle(); - } else if (edges.size() ==1) { - return edges[0]; - } else { - - CGAL_assertion(edges.size()== 2); - std::cout << "Both edges " << edges[0]->curve() << " and " - << edges[1]->curve() - << " are OK" << std::endl; - if (edges[0]->vertex()== edges[1]->opposite()->vertex()) { - throw On_vertex_exception(edges[0]->vertex()); - } else { - CGAL_assertion(edges[1]->vertex()== edges[0]->opposite()->vertex()); - throw On_vertex_exception(edges[1]->vertex()); - } - } + CGAL_assertion(0); + return Halfedge_handle(); } @@ -150,34 +112,23 @@ CGAL::Comparison_result Slice::rule_shoot_edge_vertex(const T::Sphere_point_3 & return rule_shoot_compare_SR(ep, rule, hp, p.rule_key(),p, arc_top); } else { - // rule from same sphere - CGAL_assertion(hp.is_inside()); - bool exact; - CGAL::Comparison_result debug_answer= debug_rule_shoot_answer(ep, rule, - Sds::Point(hp, hn), - exact); - if (hp.part() & hn.part() & rule.part()) { - // check rational should get this - CGAL_assertion(0); - /*CGAL::Comparison_result ret - = CGAL::compare(t_.center_c(hn.key(), rule.constant_coordinate()), - ep.simple_coordinate(rule.constant_coordinate())); - debug_rule_shoot_check(debug_answer, ret, exact);*/ - return CGAL::SMALLER; - } else if (hp.is_top()&& hn.is_top() - || hp.is_right() && hn.is_right() - /*hp.part() & hn.part() & Sds::Curve::T_BIT - || hp.part() & hn.part() & Sds::Curve::R_BIT*/) { - debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact); - return CGAL::LARGER; - } else { - debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact); - return CGAL::SMALLER; - } + // compare rational should have picked this up + CGAL_assertion(0); } } } else { - return rule_shoot_compare_SS(ep, rule, hp, p, hn); + bool tangent= (hp.key()== hn.key()); + if (tangent) { + std::pair ep3 + = t_.intersection_2_events(p.sphere_key(0), + p.sphere_key(1)); + std::cout << "Tangent." << std::endl; + CGAL::Comparison_result cr= ep3.first.compare(ep, sweep_coordinate()); + //debug_rule_shoot_check(debug_answer, cr, exact); + return cr; + } else { + return rule_shoot_compare_SS(ep, rule,p); + } } CGAL_assertion(0); @@ -187,9 +138,7 @@ CGAL::Comparison_result Slice::rule_shoot_edge_vertex(const T::Sphere_point_3 & CGAL::Comparison_result Slice::rule_shoot_compare_SS(const T::Sphere_point_3 & ep, Sds::Curve srule, - Sds::Curve arc0, - Sds::Point pt, - Sds::Curve arc1) const { + Sds::Point pt) const { bool exact; CGAL::Comparison_result debug_answer= debug_rule_shoot_answer(ep, srule, pt, @@ -197,84 +146,12 @@ CGAL::Comparison_result Slice::rule_shoot_compare_SS(const T::Sphere_point_3 & e //std::cout << "SS " << arc0 << "--" << pt << "--" << arc1 << std::endl; T::Coordinate_index C= srule.constant_coordinate(); - bool tangent= (arc0.key()==arc1.key()); - if (tangent) { - CGAL_assertion(0); - std::cout << "Tangent." << std::endl; - /*if (arc0.key() == pt.sphere_key(0)) arc1= pt.sphere(1); - else arc1= pt.sphere(0); - std::cout << "Using " << arc0 << "--" << pt << "--" << arc1 << std::endl;*/ - } - // NOTE maybe make these into index based predicates - bool hi= t_.sphere_intersects_rule(arc0.key(), srule.key(), C); - bool ohi= t_.sphere_intersects_rule(arc0.key(), srule.key(), C); - CGAL::Comparison_result hp= t_.compare_sphere_centers_c(arc0.key(), - srule.key(), C); - CGAL::Comparison_result ohp= t_.compare_sphere_centers_c(arc1.key(), - srule.key(), C); - //std::cout << "Comparisons are " << hi << " " << hp << " and " << ohi << " " << ohp << std::endl; - if (hi && !ohi || !ohi && !hi) { - debug_rule_shoot_check(debug_answer,ohp, exact); - return ohp; - } else if (ohi && !hi) { - debug_rule_shoot_check(debug_answer,hp, exact); - return hp; - } /*else if (!ohp && !hp) { - debug_rule_shoot_check(debug_answer,hp, exact); - CGAL_assertion(hp == ohp); - return ohp; - }*/ - // now they both intersect the plane - - - // NOTE hmmm, need to clean - CGAL::Bounded_side bs - = t_.bounded_side_of_sphere_on_equipower_plane_rule_line(arc0.key(), - arc1.key(), - srule.key(), - T::Coordinate_index(C), - ep); - if (bs == CGAL::ON_UNBOUNDED_SIDE){ - //std::cout << "Missed." << std::endl; - CGAL::Comparison_result ret= t_.compare_equipower_point_to_rule(arc0.key(), - arc1.key(), - srule.key(), C); - debug_rule_shoot_check(debug_answer,ret, exact); - return ret; - } else if (bs== CGAL::ON_BOUNDARY) { - //std::cout << "Equal one of them " << std::endl; - debug_rule_shoot_check(debug_answer,CGAL::EQUAL, exact); - return CGAL::EQUAL; - } else { - if (arc0.is_inside() || arc1.is_inside()) { - // arcs are at extremum of face if I am checking - CGAL::Sign sn= t_.sign_of_separating_plane_normal_c(pt.sphere_key(0), - pt.sphere_key(1), C); - CGAL_assertion(sn != CGAL::ZERO); - CGAL::Comparison_result cret= CGAL::enum_cast(sn); - debug_rule_shoot_check(debug_answer, - cret, - exact); - return cret; - } else { - //bool between =true; - CGAL::Sign dir= t_.sign_of_equipower_plane_normal_c(arc0.key(), arc1.key(), - C); - //CGAL::Comparison_result dir=CGAL::compare(eqp.orthogonal_vector()[C], 0); - int cum=0; - CGAL::Oriented_side os= t_.oriented_side_of_equipower_plane(arc0.key(), arc1.key(), ep); - if (os == CGAL::ON_POSITIVE_SIDE) ++cum; - if (dir == CGAL::POSITIVE) ++cum; - if (cum%2==0) { - debug_rule_shoot_check(debug_answer,CGAL::SMALLER, exact); - return CGAL::SMALLER; - } else { - debug_rule_shoot_check(debug_answer,CGAL::LARGER, exact); - return CGAL::LARGER; - } - } - } + CGAL::Comparison_result cr= t_.compare_sphere_sphere_at_sweep(pt.sphere_key(0), + pt.sphere_key(1), + ep, ep, C); + debug_rule_shoot_check(debug_answer, cr, exact); + return cr; } @@ -296,45 +173,32 @@ CGAL::Comparison_result Slice::rule_shoot_compare_SR(const T::Sphere_point_3 & e //CGAL_precondition(orule.is_rule()); //CGAL_precondition(srule.constant_coordinate() != orule.constant_coordinate()); - + CGAL::Comparison_result c= t_.compare_sphere_center_c(arc.key(), ep, + srule.constant_coordinate()); if (srule.is_vertical()) { - CGAL::Comparison_result c= t_.compare_sphere_centers_c(srule.key(), arc.key(), - plane_coordinate(0)); - if (c== CGAL::SMALLER && arc.is_right()) { + if (c== CGAL::LARGER && arc.is_right()) { debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact); return CGAL::LARGER; - } else if (c== CGAL::LARGER && arc.is_left()) { + } else if (c== CGAL::SMALLER && arc.is_left()) { debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact); return CGAL::SMALLER; } } else { - CGAL::Comparison_result c= t_.compare_sphere_centers_c(srule.key(), arc.key(), - plane_coordinate(1)); - if (c== CGAL::SMALLER && arc.is_top()) { - debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact); - return CGAL::LARGER; - } else if (c== CGAL::LARGER && !arc.is_top()) { + if (c== CGAL::SMALLER && arc.is_bottom()) { debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact); return CGAL::SMALLER; + } else if (c== CGAL::LARGER && arc.is_top()) { + debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact); + return CGAL::LARGER; } } - CGAL::Bounded_side intersection_side; - if (srule.constant_coordinate() == plane_coordinate(0)) { - intersection_side= t_.bounded_side_of_sphere(arc.key(), - srule.key(), - orule, - ep); - } else { - intersection_side= t_.bounded_side_of_sphere(arc.key(), - orule, - srule.key(), - ep); - } - //ep, orule, - // arc.key(), - // srule.constant_coordinate()); + CGAL::Bounded_side intersection_side + =t_.bounded_side_of_sphere_projected(arc.key(), + ep, + orule, + srule.constant_coordinate()); if (intersection_side== CGAL::ON_BOUNDARY){ @@ -418,7 +282,7 @@ void Slice::debug_rule_shoot_check(CGAL::Comparison_result check, } -bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep, +/*bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep, Sds::Curve rule, Sds::Curve a, CGAL::Comparison_result &ret) const { @@ -443,7 +307,7 @@ bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep, } } return false; -} + }*/ bool Slice::rule_shoot_compare_if_rational(const T::Sphere_point_3 & ep, @@ -452,82 +316,20 @@ bool Slice::rule_shoot_compare_if_rational(const T::Sphere_point_3 & ep, CGAL::Comparison_result &ret) const { //NT coord; if (pt.is_rule_rule()) { - ret= t_.compare_sphere_centers_c(pt.rule_key(rule.constant_coordinate()), - rule.key(), - rule.constant_coordinate()); + ret= t_.compare_sphere_center_c(pt.rule_key(rule.constant_coordinate()), + ep, + rule.constant_coordinate()); return true; } else if (pt.is_sphere_rule() && pt.rule_coordinate() == rule.constant_coordinate()) { - ret= t_.compare_sphere_centers_c(pt.rule_key(), - rule.key(), - rule.constant_coordinate()); + ret= t_.compare_sphere_center_c(pt.rule_key(), + ep, + rule.constant_coordinate()); return true; } else { - // need the vertex_handle - // actually, I don't think it was ever useful return false; } - /*if (a.is_arc() && rule_shoot_compare_if_rational_arc(ep, rule, a, ret)){ - return true; - } - if (b.is_arc() && rule_shoot_compare_if_rational_arc(ep, rule, b, ret)){ - return true; - } - - return false; - }*/ - /*std::cout << "It is rational with coordinate " - << CGAL::to_double(coord) << std::endl;*/ - - /*ret= CGAL::compare(coord, - t_.center(rule.key())[rule.constant_coordinate()]); - return true;*/ } -// compare along C where the intersection point of the two rules -// (defined by Coord C of ep) is relative to the sphere of the arc -#if 0 -CGAL::Bounded_side Slice::rule_shoot_source_side(const T::Sphere_point_3 & ep, - T::Key rule, - Sds::Curve orule, - T::Key arc_index, - int C) const { - return sc_.bounded_side_of_sphere_rule_rule_line(epi, - orule.index(), - arc_index, - C, - ep); - /*NT p[2]; - p[C]= ep.simple_coordinate(C); - p[1-C]= center_c(orule.key(), orule.constant_coordinate()); - - std::cout << "Line through " << CGAL::to_double(p[0]) - << " " << CGAL::to_double(p[1]) << std::endl; - T::Line_3 l(T::Point_3(p[0], p[1], 0), T::Vector_3(0,0,1)); - - const T::Sphere_point_3 & p0=const T::Sphere_point_3 &(sphere(arc_index), l); - CGAL::Bounded_side intersection_side; - if (!p0.is_valid()) { - intersection_side=CGAL::ON_UNBOUNDED_SIDE; - } else if (p0 == ep) { - std::cout << "Equal" << std::endl; - intersection_side= CGAL::ON_BOUNDARY; - } else { - const T::Sphere_point_3 & p1=const T::Sphere_point_3 &(sphere(arc_index), - l.opposite()); - if (p1 == ep) { - std::cout << "Equal 1" << std::endl; - intersection_side= CGAL::ON_BOUNDARY; - } - CGAL_assertion(p0 <= p1); - if (p0 < ep && p1 > ep){ - intersection_side= CGAL::ON_BOUNDED_SIDE; - } else { - intersection_side= CGAL::ON_UNBOUNDED_SIDE; - } - } - return intersection_side;*/ -} -#endif diff --git a/Arrangement_of_spheres_3/src/CGAL/makefile.internal b/Arrangement_of_spheres_3/src/CGAL/makefile.internal index 9b2c1cec3f3..b4243461fdb 100644 --- a/Arrangement_of_spheres_3/src/CGAL/makefile.internal +++ b/Arrangement_of_spheres_3/src/CGAL/makefile.internal @@ -15,7 +15,8 @@ include $(CGAL_MAKEFILE) # compiler flags #---------------------------------------------------------------------# # -fmudflapth -CXXFLAGS = -I../../include -Wall -DCGAL_USE_GMP -DCGAL_USE_GMPXX -DCGAL_USE_CORE -DCGAL_USE_QT -DQT_THREAD_SUPPORT '-I$(CGAL_INCL_DIR)' '-I$(CGAL_INCL_DIR)/CGAL/config/$(CGAL_OS_COMPILER)' '-I$(QTDIR)/include' -g -DCGAL_CHECK_EXPENSIVE -DCGAL_CHECK_EXACTNESS +CXXFLAGS = -I../../include -Wall -DCGAL_USE_GMP -DCGAL_USE_GMPXX -DCGAL_USE_CORE -DCGAL_USE_QT -DQT_THREAD_SUPPORT '-I$(CGAL_INCL_DIR)' '-I$(CGAL_INCL_DIR)/CGAL/config/$(CGAL_OS_COMPILER)' '-I$(QTDIR)/include' -g -DCGAL_CHECK_EXPENSIVE +# -DCGAL_CHECK_EXACTNESS # $(CGAL_LIB_CXXFLAGS) -I../include #