diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/diag.spheres b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/diag.spheres index 08161537906..6a320ee43b9 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/diag.spheres +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/diag.spheres @@ -1,5 +1,5 @@ -30 -6 -6 1.21 1 -30 -6.5 -5.5 1.21 1 -30 -7 -5 1.21 1 -30 -7.5 -4.5 1.21 1 -30 -8 -4 1.21 1 +30.00001 -6.0 -6.0 1.210001 1 +30.00002 -6.5 -5.5 1.210002 1 +30.00003 -7.0 -5.0 1.210003 1 +30.00004 -7.5 -4.5 1.210004 1 +30.00005 -8.0 -4.0 1.210005 1 diff --git a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/problem.spheres b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/problem.spheres index 23856d0a617..51b83eb1bdc 100644 --- a/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/problem.spheres +++ b/Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/problem.spheres @@ -1,7 +1,3 @@ -#-15.373 87.521 30.811 2.89 1 -#-16.6 87.831 31.656 2.89 1 -#-17.373 86.937 31.997 2.3104 1 --15.613 87.924 29.336 2.95 1 --16.032 89.293 29.267 2.3104 1 --16.672 87.03 28.702 2.89 1 -#2.3104 +30.00001 -6.0001 -6.0003 1.210001 1 +30.00002 -6.5002 -5.5004 1.210002 1 + diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Arrangement_of_spheres_traits_3_impl.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Arrangement_of_spheres_traits_3_impl.h index fa64256872b..7cafba96f76 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Arrangement_of_spheres_traits_3_impl.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Arrangement_of_spheres_traits_3_impl.h @@ -271,9 +271,6 @@ Arrangement_of_spheres_traits_3 CGAL_AOS3_TARG::intersection_2_events(Sphere_3_k - - - CGAL_AOS3_TEMPLATE CGAL::Comparison_result Arrangement_of_spheres_traits_3 CGAL_AOS3_TARG::compare_sphere_sphere_at_sweep(const Sphere_point_3 &t, diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section.h index 48b98203e59..4503b286314 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section.h @@ -42,7 +42,7 @@ public: } }; -struct Handle_equal{ + struct Handle_equal{ template bool operator()(H a, H b) const { return &*a == &*b; @@ -210,6 +210,11 @@ struct Handle_equal{ std::cout << " with " << *cc << std::endl; } + for (It ic= ib, fc= fb; ic != ie; ++ic, ++fc, ++cc) { + CGAL_assertion(is_valid(*ic)); + CGAL_assertion(is_valid(*fc)); + } + --num_components_; Face_handle f= (*fb)->face(); for (It ic=ib, fc= fb; ic != ie; ++ic, ++fc) { @@ -300,7 +305,9 @@ struct Handle_equal{ h=h->next(); } while (h->opposite() != *cp1); } - + for (It c= b; c !=e ; ++c){ + CGAL_assertion(is_valid(*c)); + } #endif std::cout << "Snipping out "; @@ -518,6 +525,27 @@ struct Handle_equal{ // insert the vertex so that h->opposite points to it Halfedge_handle insert_vertex(Vertex_handle vh, Halfedge_handle h); + bool is_valid(Halfedge_const_handle h) const { + for (Halfedge_const_iterator it =halfedges_begin(); it != halfedges_end(); ++it) { + if (it==h) return true; + } + return false; + } + + bool is_valid(Vertex_const_handle h) const { + for (Vertex_const_iterator it =vertices_begin(); it != vertices_end(); ++it) { + if (it==h) return true; + } + return false; + } + + bool is_valid(Face_const_handle h) const { + for (Face_const_iterator it =faces_begin(); it != faces_end(); ++it) { + if (it==h) return true; + } + return false; + } + private: diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section_impl.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section_impl.h index 808439737a9..cfe0f24135d 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section_impl.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Combinatorial_cross_section_impl.h @@ -35,6 +35,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::a_halfedge(Curve::Key k) const { CGAL_precondition(halfedges_[k.input_index()]== Halfedge_handle() || halfedges_[k.input_index()]->curve().is_inside() && halfedges_[k.input_index()]->curve().key() == k); + CGAL_precondition(halfedges_[k.input_index()]== Halfedge_handle() || is_valid(halfedges_[k.input_index()])); return halfedges_[k.input_index()]; } @@ -42,12 +43,14 @@ Combinatorial_cross_section CGAL_AOS3_TARG::a_halfedge(Curve::Key k) const { CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::set_halfedge(Halfedge_handle h){ + CGAL_precondition(h== Halfedge_handle() || is_valid(h)); halfedges_[h->curve().key().input_index()] = h; } CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::next_edge_on_rule(Halfedge_handle h) const { + CGAL_precondition(h== Halfedge_handle() || is_valid(h)); CGAL_precondition(h->curve().is_rule()); if (!h->vertex()->point().is_rule_rule()) return Halfedge_handle(); //int deg = degree(h->vertex()); @@ -66,6 +69,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::next_edge_on_rule(Halfedge_handle h) CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::next_edge_on_arc(Halfedge_handle h) const { + CGAL_precondition(h== Halfedge_handle() || is_valid(h)); CGAL_precondition(h->curve().is_arc()); Halfedge_handle r= h->next(); do { @@ -86,6 +90,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::next_edge_on_arc(Halfedge_handle h) CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::cross_edge(Halfedge_handle h) const { + CGAL_precondition(h== Halfedge_handle() || is_valid(h)); CGAL_precondition(degree(h->vertex()) == 3 || degree(h->vertex()) ==4 ); Halfedge_handle n= h->next(); if (h->curve().is_rule() && n->curve().is_rule() @@ -100,6 +105,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::cross_edge(Halfedge_handle h) const CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::is_redundant(Vertex_const_handle v) const { + CGAL_precondition(is_valid(v)); if (degree(v) != 2) return false; if (v->point().is_rule_rule()) return true; else { @@ -109,6 +115,7 @@ bool Combinatorial_cross_section CGAL_AOS3_TARG::is_redundant(Vertex_const_handl CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::is_redundant(Halfedge_const_handle v) const { + CGAL_precondition(is_valid(v)); if (v->curve().is_arc()) return false; else { return !v->vertex()->point().is_sphere_extremum() && !v->opposite()->vertex()->point().is_sphere_extremum(); @@ -117,6 +124,7 @@ bool Combinatorial_cross_section CGAL_AOS3_TARG::is_redundant(Halfedge_const_han CGAL_AOS3_TEMPLATE unsigned int Combinatorial_cross_section CGAL_AOS3_TARG::degree(Vertex_const_handle v) const { + CGAL_precondition(is_valid(v)); unsigned int r=0; Halfedge_const_handle h= v->halfedge(); do { @@ -131,7 +139,8 @@ CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::split_face(Curve c, Halfedge_handle o, Halfedge_handle d){ - + CGAL_precondition(is_valid(o)); + CGAL_precondition(is_valid(d)); CGAL_precondition(o->face() == d->face()); std::cout << "Spliting face "; write(o->face(), std::cout) << std::endl; @@ -163,6 +172,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::split_face(Curve c, CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::remove_vertex(Vertex_handle v) { + CGAL_precondition(is_valid(v)); Halfedge_handle h= v->halfedge(); std::cout << "Removing vertex " << h->vertex()->point() << " from edges " << h->curve() @@ -206,7 +216,8 @@ Combinatorial_cross_section CGAL_AOS3_TARG::remove_vertex(Vertex_handle v) { CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::move_target(Halfedge_handle h, Vertex_handle v, bool cleanup) { - + CGAL_precondition(is_valid(h)); + CGAL_precondition(is_valid(v)); v_.on_merge_faces(h); // extra work, but... std::cout << "Moving edge "; write(h, std::cout) << " to vertex "; @@ -264,6 +275,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::move_target(Halfedge_handle h, Verte CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Face_handle Combinatorial_cross_section CGAL_AOS3_TARG::join_face(Halfedge_handle h, bool clean_up) { + CGAL_precondition(is_valid(h)); std::cout << "Merging faces "; write(h->face(), std::cout) << " and "; write(h->opposite()->face(), std::cout) << std::endl; @@ -342,6 +354,8 @@ Combinatorial_cross_section CGAL_AOS3_TARG::merge_faces(Vertex_handle v) { CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::find_halfedge(Vertex_handle v, Face_handle f) { + CGAL_precondition(is_valid(v)); + CGAL_precondition(is_valid(f)); Halfedge_handle e=v->halfedge(), h=e; do { if (h->face()==f) return h; @@ -420,6 +434,9 @@ CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::insert_vertex(Vertex_handle v, Halfedge_handle h) { + CGAL_precondition(is_valid(v)); + CGAL_precondition(is_valid(h)); + Halfedge_handle i=h->opposite(); Vertex_handle vh= h->vertex(); @@ -551,6 +568,8 @@ Combinatorial_cross_section CGAL_AOS3_TARG::unpinch_bl(Halfedge_handle a, Halfed CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::connect(Halfedge_handle a, Halfedge_handle b) { + CGAL_precondition(is_valid(a)); + CGAL_precondition(is_valid(b)); std::cout << "connecting "; write(a, std::cout) << " to "; write(b, std::cout) << std::endl; @@ -562,6 +581,7 @@ CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::new_circle(Curve::Key k, Halfedge_handle vs[4]) { + //audit(true); ++num_components_; /* std::cout << "Auditing const decorator..." << std::flush; @@ -782,17 +802,20 @@ Combinatorial_cross_section CGAL_AOS3_TARG::insert_target(Curve::Key k, CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::is_in_slice(Vertex_const_handle v) const{ + CGAL_precondition(is_valid(v)); return !v->point().is_special(); } CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::is_in_slice(Halfedge_const_handle h) const{ + CGAL_precondition(is_valid(h)); return h->face() != inf_; // really could check if I am the outside edge of the box } CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::is_in_slice(Face_const_handle h) const{ + CGAL_precondition(is_valid(h)); return h != inf_ && is_in_slice(h->halfedge()); } @@ -884,6 +907,9 @@ Combinatorial_cross_section CGAL_AOS3_TARG::remove_target(Halfedge_handle ts[4], CGAL_AOS3_TEMPLATE bool Combinatorial_cross_section CGAL_AOS3_TARG::has_vertex(Face_const_handle fh, Vertex_const_handle vh) const { + CGAL_precondition(is_valid(fh)); + CGAL_precondition(is_valid(vh)); + Halfedge_const_handle h= vh->halfedge(); do { if (h->face() == fh) return true; @@ -910,6 +936,8 @@ bool Combinatorial_cross_section CGAL_AOS3_TARG::has_vertex(Face_const_handle fh }*/ CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::audit_vertex(Vertex_const_handle v, bool extra) const { + CGAL_precondition(is_valid(v)); + Point pt = v->point(); if (!is_in_slice(v)) { @@ -1193,6 +1221,7 @@ void Combinatorial_cross_section CGAL_AOS3_TARG::audit(bool extra_vertices) cons CGAL_AOS3_TEMPLATE std::ostream &Combinatorial_cross_section CGAL_AOS3_TARG::write(Vertex_const_handle h, std::ostream &out) const { + CGAL_precondition(is_valid(h)); out << h->point(); return out; } @@ -1201,6 +1230,7 @@ std::ostream &Combinatorial_cross_section CGAL_AOS3_TARG::write(Vertex_const_han CGAL_AOS3_TEMPLATE std::ostream &Combinatorial_cross_section CGAL_AOS3_TARG::write(Halfedge_const_handle h, std::ostream &out) const { + CGAL_precondition(is_valid(h)); if (h == Halfedge_const_handle()) out << "NULL"; else { out << h->opposite()->vertex()->point() << " -- " << h->curve() @@ -1213,6 +1243,7 @@ std::ostream &Combinatorial_cross_section CGAL_AOS3_TARG::write(Halfedge_const_h CGAL_AOS3_TEMPLATE std::ostream &Combinatorial_cross_section CGAL_AOS3_TARG::write(Face_const_handle f, std::ostream &out) const { + CGAL_precondition(is_valid(f)); if (f->halfedge() == Halfedge_const_handle()) out << "empty"; else { Halfedge_const_handle e= f->halfedge(), h= f->halfedge(); @@ -1458,6 +1489,7 @@ void Combinatorial_cross_section CGAL_AOS3_TARG::init_halfedge(Halfedge_handle h CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::delete_edge(Halfedge_handle h) { + CGAL_precondition(is_valid(h)); #ifndef NDEBUG if (h->curve().key().is_input()) { CGAL_assertion(halfedges_[h->curve().key().input_index()]!=h); @@ -1566,6 +1598,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::unintersect(Face_handle fn) { CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::delete_component(Vertex_handle vh) { + CGAL_precondition(is_valid(vh)); std::cout << "Deleting component from "; write(vh, std::cout) << std::endl; --num_components_; @@ -1630,6 +1663,7 @@ Combinatorial_cross_section CGAL_AOS3_TARG::delete_component(Vertex_handle vh) { CGAL_AOS3_TEMPLATE void Combinatorial_cross_section CGAL_AOS3_TARG::delete_circle(Vertex_handle v) { + CGAL_precondition(is_valid(v)); halfedges_[v->point().key().input_index()]=Halfedge_handle(); delete_component(v); } @@ -1640,6 +1674,8 @@ Combinatorial_cross_section CGAL_AOS3_TARG::delete_circle(Vertex_handle v) { CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Halfedge_handle Combinatorial_cross_section CGAL_AOS3_TARG::halfedge(Vertex_handle v, Face_handle f) const { + CGAL_precondition(is_valid(v)); + CGAL_precondition(is_valid(f)); Halfedge_handle h= v->halfedge(); do { if (h->face() == f) return h; @@ -1655,6 +1691,8 @@ CGAL_AOS3_TEMPLATE void CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::swap_curves(Halfedge_handle ha, Halfedge_handle hb) { + CGAL_precondition(is_valid(ha)); + CGAL_precondition(is_valid(hb)); Curve t= ha->curve(); Curve to= ha->opposite()->curve(); ha->set_curve(hb->curve()); @@ -1672,6 +1710,8 @@ CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Face_handle CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::intersect_arcs(Halfedge_handle hl, Halfedge_handle hk) { + CGAL_precondition(is_valid(hl)); + CGAL_precondition(is_valid(hk)); v_.on_delete_edge(hl); v_.on_delete_edge(hk); Point pkl(hk->curve(),hl->curve()); @@ -1721,7 +1761,9 @@ CGAL_AOS3_TEMPLATE CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::Face_handle CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::unintersect_arcs(Halfedge_handle hl, Halfedge_handle hk) { - + CGAL_precondition(is_valid(hl)); + CGAL_precondition(is_valid(hk)); + v_.on_delete_edge(hl); v_.on_delete_edge(hk); v_.on_delete_edge(next_edge_on_arc(hl)); @@ -1763,6 +1805,7 @@ CGAL_AOS3_TYPENAME Combinatorial_cross_section CGAL_AOS3_TARG::handle_new_edges( edges.erase(std::unique(edges.begin(), edges.end(), Handle_equal()), edges.end()); for (unsigned int i=0; i< edges.size(); ++i) { v_.on_new_edge(edges[i]); + CGAL_precondition(is_valid(edges[i])); } } diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_processor_impl.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_processor_impl.h index b0d9cb521c7..00dccb76ebd 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_processor_impl.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_processor_impl.h @@ -10,146 +10,237 @@ Event_processor CGAL_AOS3_TARG::Event_processor(CCS &cs, Traits &tr): cs_(cs), t } +CGAL_AOS3_TEMPLATE +void Event_processor CGAL_AOS3_TARG::handle_degeneracy() { + CGAL_AOS3_TYPENAME Traits::Event_point_3 ep= cs_.visitor().simulator()->current_time(); + std::cout << "Degeneracy at " << ep << std::endl; + std::vector vertices; + ICS ics(tr_, cs_); + for (CGAL_AOS3_TYPENAME CCS::Vertex_const_iterator it= cs_.vertices_begin(); + it != cs_.vertices_end(); ++it) { + if (ics.equal_point(it, ep)) { + std::cout << it->point() << std::endl; + } + } + CGAL_assertion(0); +} + + +CGAL_AOS3_TEMPLATE +void Event_processor CGAL_AOS3_TARG::check_degeneracy() { + if (cs_.visitor().simulator()->current_time() + == cs_.visitor().simulator()->next_event_time()) { + std::cout << "Degeneracy at " << cs_.visitor().simulator()->current_time() << std::endl; + throw Degeneracy(); + } +} CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::insert(Sphere_3_key k) { - ICSI icsi(tr_, cs_); try { - CGAL_AOS3_TYPENAME CCS::Face_handle f= icsi.locate(k); - //slice.new_marked_face(f); - icsi.insert(k, f); - } catch (CGAL_AOS3_TYPENAME ICSI::On_edge_exception e) { - icsi.insert(k, e.halfedge_handle()); - } catch (CGAL_AOS3_TYPENAME ICSI::On_vertex_exception v) { - icsi.insert(k, v.vertex_handle()); + check_degeneracy(); + ICSI icsi(tr_, cs_); + try { + CGAL_AOS3_TYPENAME CCS::Face_handle f= icsi.locate(k); + //slice.new_marked_face(f); + icsi.insert(k, f); + } catch (CGAL_AOS3_TYPENAME ICSI::On_edge_exception e) { + icsi.insert(k, e.halfedge_handle()); + } catch (CGAL_AOS3_TYPENAME ICSI::On_vertex_exception v) { + icsi.insert(k, v.vertex_handle()); + } + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::remove(Sphere_3_key k) { - ICSR icsr(tr_, cs_); - icsr.remove_sphere(k); + try { + check_degeneracy(); + ICSR icsr(tr_, cs_); + icsr.remove_sphere(k); + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); + } } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::intersect(Sphere_3_key k, Sphere_3_key l) { - std::vector border_faces_k, - border_faces_l, border_both; - // can I determine if I need inside or outside or both - // do point/sphere location for insertion and use that to filter further - gather_incident_faces(k, std::back_inserter(border_faces_k)); - gather_incident_faces(l, std::back_inserter(border_faces_l)); - std::set_intersection(border_faces_k.begin(), border_faces_k.end(), - border_faces_l.begin(), border_faces_l.end(), - std::back_inserter(border_both), - CGAL_AOS3_TYPENAME CCS::Handle_compare()); - CGAL_AOS3_TYPENAME CCS::Face_handle f; - if (border_both.size() > 1) { - ICS ics(tr_, cs_); + try { + check_degeneracy(); + std::vector border_faces_k, + border_faces_l, border_both; + // can I determine if I need inside or outside or both + // do point/sphere location for insertion and use that to filter further + gather_incident_faces(k, std::back_inserter(border_faces_k)); + gather_incident_faces(l, std::back_inserter(border_faces_l)); + std::set_intersection(border_faces_k.begin(), border_faces_k.end(), + border_faces_l.begin(), border_faces_l.end(), + std::back_inserter(border_both), + CGAL_AOS3_TYPENAME CCS::Handle_compare()); + CGAL_AOS3_TYPENAME CCS::Face_handle f; + if (border_both.size() > 1) { + ICS ics(tr_, cs_); - f=ics.locate(border_both.begin(), border_both.end(), - cs_.visitor().simulator()->current_time()); - //f= ics.locate(border_both.begin(), border_both.end()); - } else { - f= border_both.back(); + f=ics.locate(border_both.begin(), border_both.end(), + cs_.visitor().simulator()->current_time()); + //f= ics.locate(border_both.begin(), border_both.end()); + } else { + f= border_both.back(); + } + + // have problem with more than one vertex on face + + std::cout << "Intersecting on face "; + cs_.write(f, std::cout) << std::endl; + + Halfedge_handle hk=f->halfedge(); + do { + hk=hk->next(); + } while (!hk->curve().is_arc() || hk->curve().key() != k); + + { + Halfedge_handle hkn=hk; + do { + hkn=hkn->next(); + if (hkn == f->halfedge()) break; + if (hkn->curve().is_arc() && hkn->curve().key() ==k) { + Degeneracy d; + d.new_edge(hkn); + d.new_edge(hk); + throw d; + } + } while (true); + } + + Halfedge_handle hl=f->halfedge(); + do { + hl=hl->next(); + } while (!hl->curve().is_arc() || hl->curve().key() != l); + + { + Halfedge_handle hln=hl; + do { + hln=hln->next(); + if (hln == f->halfedge()) break; + if (hln->curve().is_arc() && hln->curve().key() == l) { + Degeneracy d; + d.new_edge(hln); + d.new_edge(hl); + throw d; + } + } while (true); + } + + std::cout << "Intersection edges are "; + cs_.write(hk, std::cout) << " and "; + cs_.write(hl, std::cout) << std::endl; + + cs_.intersect_arcs(hl, hk); + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } - - std::cout << "Intersecting on face "; - cs_.write(f, std::cout) << std::endl; - - Halfedge_handle hk=f->halfedge(); - do { - hk=hk->next(); - } while (!hk->curve().is_arc() || hk->curve().key() != k); - Halfedge_handle hl=f->halfedge(); - do { - hl=hl->next(); - } while (!hl->curve().is_arc() || hl->curve().key() != l); - - std::cout << "Intersection edges are "; - cs_.write(hk, std::cout) << " and "; - cs_.write(hl, std::cout) << std::endl; - - cs_.intersect_arcs(hl, hk); } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::unintersect(Sphere_3_key k, Sphere_3_key l) { - CGAL_AOS3_TYPENAME CCS::Halfedge_handle hl, hk= cs_.a_halfedge(k)->opposite(); - CGAL_AOS3_TYPENAME CCS::Halfedge_handle khs=hk; - do { - if (hk->next()->curve().key() == l - && hk->next()->next() == hk){ - hl= hk->next(); - break; - } else if (hk->opposite()->next()->curve().key() ==l - && hk->opposite()->next()->next() == hk->opposite()) { - hl= hk->opposite()->next(); - hk= hk->opposite(); - break; + try { + check_degeneracy(); + CGAL_AOS3_TYPENAME CCS::Halfedge_handle hl, hk= cs_.a_halfedge(k)->opposite(); + CGAL_AOS3_TYPENAME CCS::Halfedge_handle khs=hk; + do { + if (hk->next()->curve().key() == l + && hk->next()->next() == hk){ + hl= hk->next(); + break; + } else if (hk->opposite()->next()->curve().key() ==l + && hk->opposite()->next()->next() == hk->opposite()) { + hl= hk->opposite()->next(); + hk= hk->opposite(); + break; + } + hk= cs_.next_edge_on_arc(hk); + } while (hk != khs); + if (hl == Halfedge_handle() || hk== Halfedge_handle()) { + throw Degeneracy(); } - hk= cs_.next_edge_on_arc(hk); - } while (hk != khs); - // find shared face, f, and edges hk and hl - cs_.unintersect_arcs(hl, hk); + // find shared face, f, and edges hk and hl + cs_.unintersect_arcs(hl, hk); + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); + } } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::intersect(Sphere_3_key k, Sphere_3_key l, Sphere_3_key m) { - Halfedge_handle hl, hm, hk= cs_.a_halfedge(k)->opposite(); - Halfedge_handle khs=hk; - do { - if (hk->next()->curve().key() == l - && hk->next()->next()->curve().key() == m - && hk->next()->next()->next() == hk){ - hl= hk->next(); - hm= hk->next()->next(); - break; - } else if (hk->next()->curve().key() == m - && hk->next()->next()->curve().key() == l - && hk->next()->next()->next() == hk){ - hl= hk->next()->next(); - hm= hk->next(); - break; - } else if (hk->opposite()->next()->curve().key() == l - && hk->opposite()->next()->next()->curve().key() == m - && hk->opposite()->next()->next()->next() == hk){ - hk=hk->opposite(); - hl= hk->next(); - hm= hk->next()->next(); - break; - } else if (hk->opposite()->next()->curve().key() == m - && hk->opposite()->next()->next()->curve().key() == l - && hk->opposite()->next()->next()->next() == hk){ - hk=hk->opposite(); - hl= hk->next()->next(); - hm= hk->next(); - break; - } - hk= cs_.next_edge_on_arc(hk); - } while (hk != khs); + try { + check_degeneracy(); + Halfedge_handle hl, hm, hk= cs_.a_halfedge(k)->opposite(); + Halfedge_handle khs=hk; + do { + if (hk->next()->curve().key() == l + && hk->next()->next()->curve().key() == m + && hk->next()->next()->next() == hk){ + hl= hk->next(); + hm= hk->next()->next(); + break; + } else if (hk->next()->curve().key() == m + && hk->next()->next()->curve().key() == l + && hk->next()->next()->next() == hk){ + hl= hk->next()->next(); + hm= hk->next(); + break; + } else if (hk->opposite()->next()->curve().key() == l + && hk->opposite()->next()->next()->curve().key() == m + && hk->opposite()->next()->next()->next() == hk){ + hk=hk->opposite(); + hl= hk->next(); + hm= hk->next()->next(); + break; + } else if (hk->opposite()->next()->curve().key() == m + && hk->opposite()->next()->next()->curve().key() == l + && hk->opposite()->next()->next()->next() == hk){ + hk=hk->opposite(); + hl= hk->next()->next(); + hm= hk->next(); + break; + } + hk= cs_.next_edge_on_arc(hk); + } while (hk != khs); + if (hk == Halfedge_handle()) { + throw Degeneracy(); + } - Halfedge_handle hn= cs_.next_edge_on_arc(hk)->opposite(); - Halfedge_handle hp= cs_.next_edge_on_arc(hk); + Halfedge_handle hn= cs_.next_edge_on_arc(hk)->opposite(); + Halfedge_handle hp= cs_.next_edge_on_arc(hk); - Halfedge_handle hnt= hn->next()->next(); - Halfedge_handle hpt= hp->prev()->prev(); + Halfedge_handle hnt= hn->next()->next(); + Halfedge_handle hpt= hp->prev()->prev(); - Vertex_handle nvhn= cs_.new_vertex(Point(hk->curve(), hnt->curve())); - Vertex_handle nvhp= cs_.new_vertex(Point(hk->curve(), hpt->curve())); + Vertex_handle nvhn= cs_.new_vertex(Point(hk->curve(), hnt->curve())); + Vertex_handle nvhp= cs_.new_vertex(Point(hk->curve(), hpt->curve())); - cs_.insert_vertex(nvhn, hnt); - cs_.insert_vertex(nvhp, hpt); - cs_.move_target(hp, nvhp, true); - cs_.move_target(hn, nvhn, true); + cs_.insert_vertex(nvhn, hnt); + cs_.insert_vertex(nvhp, hpt); + cs_.move_target(hp, nvhp, true); + cs_.move_target(hn, nvhn, true); - cs_.split_face(hk->curve(), hn->next()->opposite(), - hp->prev()->opposite()->prev()); - cs_.join_face(hk, true); + cs_.split_face(hk->curve(), hn->next()->opposite(), + hp->prev()->opposite()->prev()); + cs_.join_face(hk, true); + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); + } } CGAL_AOS3_TEMPLATE @@ -160,52 +251,64 @@ void Event_processor CGAL_AOS3_TARG::unintersect(Sphere_3_key k, Sphere_3_key l, CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::process_aar(Halfedge_handle h) { - clear_event(h); - CGAL_precondition(h->vertex()->point().is_sphere_rule()); - Halfedge_handle rule= cs_.cross_edge(h); - CGAL_assertion(rule->vertex() ==h->vertex()); - Halfedge_handle target; - if (h->next() == rule->opposite()) { - target= h->prev(); - } else { - target=h->opposite()->next(); + try { + check_degeneracy(); + clear_event(h); + CGAL_precondition(h->vertex()->point().is_sphere_rule()); + Halfedge_handle rule= cs_.cross_edge(h); + CGAL_assertion(rule->vertex() ==h->vertex()); + Halfedge_handle target; + if (h->next() == rule->opposite()) { + target= h->prev(); + } else { + target=h->opposite()->next(); + } + Point pt= h->vertex()->point(); + Vertex_handle nv= cs_.new_vertex(pt); + cs_.insert_vertex(nv, target); + cs_.move_target(rule, nv, true); + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } - Point pt= h->vertex()->point(); - Vertex_handle nv= cs_.new_vertex(pt); - cs_.insert_vertex(nv, target); - cs_.move_target(rule, nv, true); } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::process_rar(Halfedge_handle h) { - clear_event(h); - if (h->opposite()->vertex()->point().is_sphere_extremum()) { - h=h->opposite(); - } - CGAL_assertion(cs_.degree(h->vertex())==3); - CGAL_assertion(cs_.degree(h->opposite()->vertex())==3); - // pick which to move, make sure don't move extremum - Halfedge_handle move_rule= cs_.cross_edge(h->opposite())->opposite(); - CGAL_assertion(move_rule->vertex() == h->opposite()->vertex()); - Halfedge_handle onto_edge= cs_.cross_edge(h); - if (onto_edge->face() != move_rule->face() - && onto_edge->face() != move_rule->opposite()->face()){ - onto_edge= onto_edge->opposite(); - } - Halfedge_handle other_rule= cs_.cross_edge(h); - // handle them not being on the same side - if (onto_edge->face() != move_rule->face() - && onto_edge->face() != move_rule->opposite()->face()){ - Vertex_handle ov= move_rule->vertex(); - cs_.move_target(move_rule, other_rule->vertex(), false); - cs_.move_target(other_rule, ov, true); - } else { - Point npt(move_rule->curve(), other_rule->curve()); - Vertex_handle nv= cs_.new_vertex(npt); - cs_.insert_vertex(nv, other_rule); - cs_.move_target(move_rule, nv, true); + try { + check_degeneracy(); + clear_event(h); + if (h->opposite()->vertex()->point().is_sphere_extremum()) { + h=h->opposite(); + } + CGAL_assertion(cs_.degree(h->vertex())==3); + CGAL_assertion(cs_.degree(h->opposite()->vertex())==3); + // pick which to move, make sure don't move extremum + Halfedge_handle move_rule= cs_.cross_edge(h->opposite())->opposite(); + CGAL_assertion(move_rule->vertex() == h->opposite()->vertex()); + Halfedge_handle onto_edge= cs_.cross_edge(h); + if (onto_edge->face() != move_rule->face() + && onto_edge->face() != move_rule->opposite()->face()){ + onto_edge= onto_edge->opposite(); + } + Halfedge_handle other_rule= cs_.cross_edge(h); + // handle them not being on the same side + if (onto_edge->face() != move_rule->face() + && onto_edge->face() != move_rule->opposite()->face()){ + Vertex_handle ov= move_rule->vertex(); + cs_.move_target(move_rule, other_rule->vertex(), false); + cs_.move_target(other_rule, ov, true); + } else { + Point npt(move_rule->curve(), other_rule->curve()); + Vertex_handle nv= cs_.new_vertex(npt); + cs_.insert_vertex(nv, other_rule); + cs_.move_target(move_rule, nv, true); + } + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } } @@ -219,30 +322,32 @@ void Event_processor CGAL_AOS3_TARG::clear_event(Halfedge_handle h) { CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::process_arr(Halfedge_handle h) { - clear_event(h); - if (!h->vertex()->point().is_rule_rule()) h= h->opposite(); - CGAL_assertion(h->vertex()->point().is_rule_rule()); - Halfedge_handle xr= cs_.cross_edge(h)->opposite(); // inward pointing - //CGAL_assertion(xr->vertex() == h->vertex()); - std::cout << "XR is "; - cs_.write(xr, std::cout) << std::endl; + try { + check_degeneracy(); + clear_event(h); + if (!h->vertex()->point().is_rule_rule()) h= h->opposite(); + CGAL_assertion(h->vertex()->point().is_rule_rule()); + Halfedge_handle xr= cs_.cross_edge(h)->opposite(); // inward pointing + //CGAL_assertion(xr->vertex() == h->vertex()); + std::cout << "XR is "; + cs_.write(xr, std::cout) << std::endl; - ICS ics(tr_, cs_); - ics.roll_back_rule(cs_.visitor().simulator()->current_time(), xr); - Halfedge_handle nxr= cs_.next_edge_on_rule(xr); - if (nxr != Halfedge_handle()) { - if (cs_.is_redundant(xr)) { - cs_.join_face(xr, true); - xr= nxr->opposite(); - std::cout << "XR is "; - cs_.write(xr, std::cout) << std::endl; - nxr=Halfedge_handle(); - } else { - nxr= nxr->opposite(); - std::cout << "NXR is "; - cs_.write(nxr, std::cout) << std::endl; + ICS ics(tr_, cs_); + ics.roll_back_rule(cs_.visitor().simulator()->current_time(), xr); + Halfedge_handle nxr= cs_.next_edge_on_rule(xr); + if (nxr != Halfedge_handle()) { + if (cs_.is_redundant(xr)) { + cs_.join_face(xr, true); + xr= nxr->opposite(); + std::cout << "XR is "; + cs_.write(xr, std::cout) << std::endl; + nxr=Halfedge_handle(); + } else { + nxr= nxr->opposite(); + std::cout << "NXR is "; + cs_.write(nxr, std::cout) << std::endl; + } } - } else { if (cs_.is_redundant(xr)) { cs_.join_face(xr, true); } else { @@ -262,133 +367,148 @@ void Event_processor CGAL_AOS3_TARG::process_arr(Halfedge_handle h) { cs_.move_target(xr, nv, true); } while ((xr=nxr) != Halfedge_handle()); // do other side if needed } + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } - // could fail if there is a degeneracy - // now the cross rule should be a t pointing on direction; + // now the cross rule should be a t pointing on direction; } CGAL_AOS3_TEMPLATE void Event_processor CGAL_AOS3_TARG::process_aae(Sphere_3_key k, Sphere_3_key l, int i) { - int index=i+1; - if (tr_.compare_sphere_center_c(k, cs_.visitor().simulator()->current_time(), plane_coordinate(1-i))== CGAL::LARGER) { - index+=2; - } - std::cout << "Extremum event between " << k << " and " << l << " on rules "; - if (i==0) std::cout << "TB"; - else std::cout << "LR"; - std::cout << " " << index << std::endl; + try { + check_degeneracy(); + int index=i+1; + if (tr_.compare_sphere_center_c(k, cs_.visitor().simulator()->current_time(), plane_coordinate(1-i))== CGAL::LARGER) { + index+=2; + } + std::cout << "Extremum event between " << k << " and " << l << " on rules "; + if (i==0) std::cout << "TB"; + else std::cout << "LR"; + std::cout << " " << index << std::endl; - // find outwards pointing rule edge - Halfedge_handle hk= cs_.a_halfedge(k); - while (!hk->vertex()->point().is_sphere_extremum() - || hk->vertex()->point().is_sphere_extremum() - && hk->vertex()->point().sphere_extremum_index().index() != index) { - cs_.write(hk, std::cout) << " " << hk->vertex()->point().is_sphere_extremum() - << " " << (hk->vertex()->point().is_sphere_extremum()? hk->vertex()->point().sphere_extremum_index().index(): -1) << std::endl; - hk= cs_.next_edge_on_arc(hk); - } - Halfedge_handle rule= hk->opposite()->prev()->opposite(); // outward rule + // find outwards pointing rule edge + Halfedge_handle hk= cs_.a_halfedge(k); + while (!hk->vertex()->point().is_sphere_extremum() + || hk->vertex()->point().is_sphere_extremum() + && hk->vertex()->point().sphere_extremum_index().index() != index) { + cs_.write(hk, std::cout) << " " << hk->vertex()->point().is_sphere_extremum() + << " " << (hk->vertex()->point().is_sphere_extremum()? + hk->vertex()->point().sphere_extremum_index().index() + : -1) << std::endl; + hk= cs_.next_edge_on_arc(hk); + } + Halfedge_handle rule= hk->opposite()->prev()->opposite(); // outward rule - std::cout << index << " hk is "; - cs_.write(hk, std::cout) << " and rule is "; - cs_.write(rule, std::cout) << std::endl; + std::cout << index << " hk is "; + cs_.write(hk, std::cout) << " and rule is "; + cs_.write(rule, std::cout) << std::endl; - // check if it ends on arc of l - if (rule->vertex()->point().is_sphere_rule() - && rule->vertex()->point().sphere_key() == l) { + // check if it ends on arc of l + if (rule->vertex()->point().is_sphere_rule() + && rule->vertex()->point().sphere_key() == l) { - Halfedge_handle new_extremum; // inside arc - Curve base_curve, rule_curve=rule->curve(); - Point pt= hk->vertex()->point(); - if (rule->opposite()->prev()->prev() != hk->opposite()) { - new_extremum=cs_.next_edge_on_arc(cs_.next_edge_on_arc(hk)); - base_curve= hk->curve(); - } else { - CGAL_assertion(rule->opposite()->prev()->prev() == hk->opposite()); - base_curve= hk->next()->curve(); - new_extremum=cs_.next_edge_on_arc(hk->opposite())->opposite(); - } + Halfedge_handle new_extremum; // inside arc + Curve base_curve, rule_curve=rule->curve(); + Point pt= hk->vertex()->point(); + if (rule->opposite()->prev()->prev() == hk->opposite()) { + base_curve= hk->next()->curve(); + new_extremum=cs_.next_edge_on_arc(hk->opposite())->opposite(); + } else { + if (rule->opposite()->prev()->prev() != hk->opposite()) { + throw Degeneracy(); + } + new_extremum=cs_.next_edge_on_arc(cs_.next_edge_on_arc(hk)); + base_curve= hk->curve(); + } - std::cout << "New extremum is "; - cs_.write(new_extremum, std::cout) << " and base curve is " << base_curve << std::endl; + std::cout << "New extremum is "; + cs_.write(new_extremum, std::cout) << " and base curve is " << base_curve << std::endl; - CGAL_postcondition(new_extremum->curve().is_inside()); - cs_.set_curve(hk, base_curve); - cs_.set_curve(cs_.next_edge_on_arc(hk), base_curve); + CGAL_postcondition(new_extremum->curve().is_inside()); + cs_.set_curve(hk, base_curve); + cs_.set_curve(cs_.next_edge_on_arc(hk), base_curve); - cs_.join_face(rule, true); - //cs_.remove_vertex(v0, ->curve()); - //cs_.remove_vertex(v1); - Vertex_handle nvh= cs_.new_vertex(pt); - Halfedge_handle nh=cs_.insert_vertex(nvh, new_extremum); - cs_.set_curve(nh->opposite(), cs_.next_edge_on_arc(nh->opposite())->curve()); - cs_.set_curve(nh->next(), cs_.next_edge_on_arc(nh->next())->curve()); - - ICS ics(tr_, cs_); - Face_handle f; - if (nh->curve().is_inside()) { - nh=nh->opposite()->prev(); - } - f= nh->face(); - Halfedge_handle hd= ics.find_rule_vertex(cs_.visitor().simulator()->current_time(), - f, - rule_curve); - cs_.split_face(rule_curve, nh, hd); - } else { - ICS ics(tr_, cs_); - ics.roll_back_rule(cs_.visitor().simulator()->current_time(), rule); - - Halfedge_handle new_extremum, new_target; // inside arc - Curve base_curve, rule_curve=rule->curve(); - Point pt= hk->vertex()->point(); - if (rule->prev()->prev()->curve().is_arc() - && rule->prev()->prev()->curve().key() == l) { - new_extremum=cs_.next_edge_on_arc(cs_.next_edge_on_arc(hk)); - base_curve= hk->curve(); - new_target= rule->prev()->prev()->opposite(); - } else { - base_curve= cs_.next_edge_on_arc(hk)->curve(); - new_extremum=cs_.next_edge_on_arc(hk->opposite())->opposite(); - new_target= rule->opposite()->next()->next()->opposite(); - } - - std::cout << "New extremum is "; - cs_.write(new_extremum, std::cout) << " and base curve is " << base_curve - << " and new target is "; - cs_.write(new_target, std::cout) << std::endl; - - - CGAL_assertion(new_target->curve().key() ==l); - CGAL_postcondition(new_extremum->curve().is_inside()); - cs_.set_curve(hk, base_curve); - cs_.set_curve(cs_.next_edge_on_arc(hk), base_curve); - - cs_.join_face(rule, true); - //cs_.remove_vertex(v0, ->curve()); - //cs_.remove_vertex(v1); - Halfedge_handle nh, nth; - { + cs_.join_face(rule, true); + //cs_.remove_vertex(v0, ->curve()); + //cs_.remove_vertex(v1); Vertex_handle nvh= cs_.new_vertex(pt); - nh=cs_.insert_vertex(nvh, new_extremum); + Halfedge_handle nh=cs_.insert_vertex(nvh, new_extremum); cs_.set_curve(nh->opposite(), cs_.next_edge_on_arc(nh->opposite())->curve()); cs_.set_curve(nh->next(), cs_.next_edge_on_arc(nh->next())->curve()); - if (nh->curve().is_inside()) nh=nh->opposite(); + + ICS ics(tr_, cs_); + Face_handle f; + if (nh->curve().is_inside()) { + nh=nh->opposite()->prev(); + } + f= nh->face(); + Halfedge_handle hd= ics.find_rule_vertex(cs_.visitor().simulator()->current_time(), + f, + rule_curve); + cs_.split_face(rule_curve, nh, hd); + } else { + ICS ics(tr_, cs_); + ics.roll_back_rule(cs_.visitor().simulator()->current_time(), rule); + + Halfedge_handle new_extremum, new_target; // inside arc + Curve base_curve, rule_curve=rule->curve(); + Point pt= hk->vertex()->point(); + if (rule->prev()->prev()->curve().is_arc() + && rule->prev()->prev()->curve().key() == l) { + new_extremum=cs_.next_edge_on_arc(cs_.next_edge_on_arc(hk)); + base_curve= hk->curve(); + new_target= rule->prev()->prev()->opposite(); + } else { + if (!rule->opposite()->next()->next()->curve().is_arc() + || rule->opposite()->next()->next()->curve().key() != l) { + throw Degeneracy(); + } + base_curve= cs_.next_edge_on_arc(hk)->curve(); + new_extremum=cs_.next_edge_on_arc(hk->opposite())->opposite(); + new_target= rule->opposite()->next()->next()->opposite(); + } + + std::cout << "New extremum is "; + cs_.write(new_extremum, std::cout) << " and base curve is " << base_curve + << " and new target is "; + cs_.write(new_target, std::cout) << std::endl; + + + CGAL_assertion(new_target->curve().key() ==l); + CGAL_postcondition(new_extremum->curve().is_inside()); + cs_.set_curve(hk, base_curve); + cs_.set_curve(cs_.next_edge_on_arc(hk), base_curve); + + cs_.join_face(rule, true); + //cs_.remove_vertex(v0, ->curve()); + //cs_.remove_vertex(v1); + Halfedge_handle nh, nth; + { + Vertex_handle nvh= cs_.new_vertex(pt); + nh=cs_.insert_vertex(nvh, new_extremum); + cs_.set_curve(nh->opposite(), cs_.next_edge_on_arc(nh->opposite())->curve()); + cs_.set_curve(nh->next(), cs_.next_edge_on_arc(nh->next())->curve()); + if (nh->curve().is_inside()) nh=nh->opposite(); + } + { + Vertex_handle ntvh= cs_.new_vertex(Point(rule_curve, new_target->curve())); + nth=cs_.insert_vertex(ntvh, new_extremum); + if (nth->curve().is_inside()) nth=nth->opposite(); + } + cs_.split_face(rule_curve, nh, nth); } - { - Vertex_handle ntvh= cs_.new_vertex(Point(rule_curve, new_target->curve())); - nth=cs_.insert_vertex(ntvh, new_extremum); - if (nth->curve().is_inside()) nth=nth->opposite(); - } - cs_.split_face(rule_curve, nh, nth); + // if so then remove, edge and both vertices, check edge label, add edge on other side, shoot rule + // else roll back rule, remove rule/vertex, add rule on other side + } catch (Degeneracy) { + std::cout << "Degeneracy" << std::endl; + handle_degeneracy(); } - // if so then remove, edge and both vertices, check edge label, add edge on other side, shoot rule - // else roll back rule, remove rule/vertex, add rule on other side - } CGAL_AOS3_END_INTERNAL_NAMESPACE diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_visitor_impl.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_visitor_impl.h index 3eb9088d12e..fcb6f5a945a 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_visitor_impl.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Event_visitor_impl.h @@ -18,11 +18,11 @@ void Event_visitor CGAL_AOS3_TARG::initialize() { Event_key k; if (ep.first >= sim_->current_time()) { k= sim_->new_event(ep.first, - CGAL_AOS3_TYPENAME Event_processor::Insert_event(j_, *it)); + CGAL_AOS3_TYPENAME Event_processor::I_event(j_, *it)); } if (ep.second >= sim_->current_time()) { k= sim_->new_event(ep.second, - CGAL_AOS3_TYPENAME Event_processor::Remove_event(j_, *it)); + CGAL_AOS3_TYPENAME Event_processor::R_event(j_, *it)); } if (k != Event_key() && k != sim_->null_event()) { free_events_.insert(k); @@ -297,13 +297,13 @@ void Event_visitor CGAL_AOS3_TARG::process_pair(Sphere_3_key a, CGAL_assertion(ep.first <= ep.second); if (ep.first >= sim_->current_time()) { - Event_key k= sim_->new_event(ep.first, CGAL_AOS3_TYPENAME Event_processor::Intersect_event(j_, a,b)); + Event_key k= sim_->new_event(ep.first, CGAL_AOS3_TYPENAME Event_processor::I2_event(j_, a,b)); if (k != Event_key() && k != sim_->null_event()) { free_events_.insert(k); } } if (ep.second >= sim_->current_time()) { - Event_key k =sim_->new_event(ep.second, CGAL_AOS3_TYPENAME Event_processor::Unintersect_event(j_, a,b)); + Event_key k =sim_->new_event(ep.second, CGAL_AOS3_TYPENAME Event_processor::U2_event(j_, a,b)); if (k != Event_key() && k != sim_->null_event()) { free_events_.insert(k); } @@ -357,13 +357,13 @@ void Event_visitor CGAL_AOS3_TARG::process_triple(Halfedge_handle h) { if (ep.first.is_valid()) { CGAL_assertion(ep.first <= ep.second); if (ep.first >= sim_->current_time()) { - Event_key k= sim_->new_event(ep.first, CGAL_AOS3_TYPENAME Event_processor::Intersect_3_event(j_, a,b,c)); + Event_key k= sim_->new_event(ep.first, CGAL_AOS3_TYPENAME Event_processor::I3_event(j_, a,b,c)); if (k != Event_key() && k != sim_->null_event()) { free_events_.insert(k); } } if (ep.second >= sim_->current_time()) { - Event_key k= sim_->new_event(ep.second, CGAL_AOS3_TYPENAME Event_processor::Unintersect_3_event(j_, a,b,c)); + Event_key k= sim_->new_event(ep.second, CGAL_AOS3_TYPENAME Event_processor::U3_event(j_, a,b,c)); if (k != Event_key() && k != sim_->null_event()) { free_events_.insert(k); } diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section.h index 1f9b40ad91c..996f8bea256 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section.h @@ -137,7 +137,10 @@ public: cs_.write(faces[0], std::cerr) << std::endl; cs_.write(faces[1], std::cerr) << std::endl; - CGAL_assertion(0); + CGAL_AOS3_TYPENAME Traits::Degeneracy_exception d; + d.new_face(faces[0]); + d.new_face(faces[1]); + throw d; } else { CGAL_assertion(!faces.empty()); CGAL_AOS3_TYPENAME CS::Halfedge_handle h= faces[0]->halfedge(); @@ -158,8 +161,14 @@ public: for (unsigned int i=0; i< faces.size(); ++i){ cs_.write(faces[i], std::cerr ) << std::endl; } - CGAL_assertion(0); - + { + CGAL_AOS3_TYPENAME Traits::Degeneracy_exception d; + for (unsigned int i=0; i< faces.size(); ++i){ + d.new_face(faces[i]); + } + + throw d; + } } CGAL_assertion(0); @@ -266,10 +275,43 @@ public: roll_back_rule(const CGAL_AOS3_TYPENAME Traits::Sphere_point_3 &t, CGAL_AOS3_TYPENAME CS::Halfedge_handle cur); + bool + equal_points(CGAL_AOS3_TYPENAME CS::Vertex_handle vh, + const CGAL_AOS3_TYPENAME Traits::Sphere_point_3 &pt) const { + if (vh->point().is_rule_rule()) { + for (unsigned int i=0; i< 2; ++i) { + Coordinate_index ci=plane_coordinate(i); + if (pt.compare(tr_.sphere_3(vh->point().rule_key(ci)).center()[ci.index()], + ci) != CGAL::EQUAL) return false; + } + return true; + } else if (vh->point().is_sphere_rule()) { + Coordinate_index ci=vh->point().rule_coordinate(); + if (pt.compare(tr_.sphere_3(vh->point().rule_key()).center()[ci.index()], + ci) != CGAL::EQUAL) return false; + // hmmmm, now what + if (tr_.bounded_side_of_sphere_projected(pt, vh->point().sphere_key(), + vh->point().rule_key(), + pt, + vh->point().rule_coordinate()) != CGAL::ON_BOUNDARY) + return false; + return true; + } else { + if ( tr_.compare_sphere_sphere_at_sweep(pt, vh->point().sphere_key(0) , + vh->point().sphere_key(1), pt, plane_coordinate(0)) != CGAL::EQUAL) + return false; + if ( tr_.compare_sphere_sphere_at_sweep(pt, vh->point().sphere_key(0), + vh->point().sphere_key(1), pt, plane_coordinate(1)) != CGAL::EQUAL) + return false; + return true; + } + } protected: + + bool locate_point_check_face(const CGAL_AOS3_TYPENAME Traits::Event_point_3 &k, CGAL_AOS3_TYPENAME CS::Face_const_handle it, std::vector &locations) const ; diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_impl.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_impl.h index b59c3e446eb..a81cf24866c 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_impl.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_impl.h @@ -576,6 +576,10 @@ Irrational_cross_section CGAL_AOS3_TARG ::roll_back_rule(const CGAL_AOS3_TYPENAM bits.push_back(cs_.next_edge_on_rule(bits.back())); } while (bits.back() != CGAL_AOS3_TYPENAME CS::Halfedge_handle()); bits.pop_back(); + for (unsigned int i=0; i< bits.size(); ++i) { + cs_.write(bits[i], std::cout) << ", "; + } + std::cout << std::endl; if (bits.size() > 1) { diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_removal.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_removal.h index ad7c86fb529..a21e4cfd516 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_removal.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Irrational_cross_section_removal.h @@ -45,7 +45,9 @@ public: P::cs_.write(vertices[i], std::cout) << " is vertex " << i << std::endl; h= h->next(); } while (h != P::cs_.a_halfedge(k)); - CGAL_assertion(deg==4); + if (deg != 4) { + throw CGAL_AOS3_TYPENAME Traits::Degeneracy_exception(); + } } // roll in each until I have a target in a face diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_3_table.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_3_table.h index 6a745d82a25..87fb9aead75 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_3_table.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_3_table.h @@ -133,6 +133,37 @@ public: void initialize_1(); void initialize_2(); + + /* + event data, not really appropriate, but it is the best place to put it + */ + struct Event_pair_data { + typedef CGAL_AOS3_TYPENAME Traits::Event_point_3 EP; + Event_pair_data():i_(-1){} + int index_; + EP events_[2]; + }; + struct Triple_data { + Event_pair_data srr_events_; + }; + struct UPair_data { + typedef CGAL_AOS3_TYPENAME Traits::Event_point_3 EP; + struct KC_pair{ + typedef KC_pair This; + KC_pair(Key k, Coordinate_index c): k_(k), c_(c){}; + Key k_; + Coordinate_index c_; + CGAL_COMPARISONS_2(k_, c_); + }; + + std::map cxr_events_; + }; + + + + std::map triple_data_; + std::map upair_data_; + Spheres spheres_; bool has_temp_; diff --git a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_key.h b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_key.h index dd89237ffd8..12b5f00775a 100644 --- a/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_key.h +++ b/Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_key.h @@ -32,9 +32,56 @@ struct Sphere_key{ static Sphere_key bl_key() {return Sphere_key(BL);} static Sphere_key tr_key() {return Sphere_key(TR);} + struct UPair{ + typedef UPair This; + UPair(Sphere_3_key a, + Sphere_3_key b): a_(std::min(a,b)), + b_(std::max(a,b)){} + CGAL_COMPARISONS_2(a_,b_); + Sphere_3_key a_,b_; + }; + + struct UTriple{ + typedef UTriple This; + UTriple(Sphere_3_key a, + Sphere_3_key b, + Sphere_3_key c): a_(std::min(a,std::min(b,c))), + b_(std::max(a,std::min(b,c))), + c_(std::max(a,std::max(b,c))){ + CGAL_assertion(a_!= b_ && b_ != c_ && c_ != a_); + } + CGAL_COMPARISONS_3(a_, b_, c_); + + Sphere_3_key a_,b_,c_; + }; + + struct Pair{ + typedef UPair This; + Pair(Sphere_3_key a, + Sphere_3_key b): a_(a), + b_(b){} + CGAL_COMPARISONS_2(a_,b_); + Sphere_3_key a_,b_; + }; + + struct Triple{ + typedef Triple This; + Triple(Sphere_3_key a, + Sphere_3_key b, + Sphere_3_key c): a_(a), b_(b), c_(c)){ + CGAL_assertion(a_!= b_ && b_ != c_ && c_ != a_); + } + CGAL_COMPARISONS_3(a_, b_, c_); + + Sphere_3_key a_,b_,c_; + }; + int id_; }; + + + CGAL_OUTPUT(Sphere_key); CGAL_AOS3_END_INTERNAL_NAMESPACE #endif 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 367b63b87ca..4b0db7f4c47 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 @@ -27,6 +27,15 @@ CGAL_BEGIN_NAMESPACE template #endif struct Arrangement_of_spheres_traits_3 { + struct Degeneracy_exception { + + template + void new_face(T){} + + template + void new_edge(T){} + }; + #ifdef CGAL_AOS3_USE_TEMPLATES typedef Arrangement_of_spheres_traits_3 This; typedef GT Geom_traits; @@ -136,10 +145,16 @@ struct Arrangement_of_spheres_traits_3 { Event_pair sphere_intersect_rule_events(Sphere_3_key a, Sphere_3_key r, Coordinate_index C) const; */ - Event_pair circle_cross_rule_events(Sphere_3_key a, Sphere_3_key b, - Sphere_3_key rs, Coordinate_index C) const; + + void advance_circle_cross_rule_event(Sphere_3_key a, Sphere_3_key b, + Sphere_3_key rs, Coordinate_index C); - Event_pair sphere_intersect_rule_rule_events(Sphere_3_key s, Sphere_3_key rx, Sphere_3_key ry) const; + Event_point_3 circle_cross_rule_event(Sphere_3_key a, Sphere_3_key b, + Sphere_3_key rs, Coordinate_index C) const; + + void advance_sphere_intersect_rule_rule_event(Sphere_3_key s, Sphere_3_key rx, Sphere_3_key ry); + + Event_point_3 sphere_intersect_rule_rule_event(Sphere_3_key s, Sphere_3_key rx, Sphere_3_key ry) const; // not really used /*Quadratic_NT intersection_c(Sphere_3_key s, Line_3 l, Coordinate_index C) const;*/ @@ -206,6 +221,8 @@ struct Arrangement_of_spheres_traits_3 { const Sphere_point_3 &pt, Coordinate_index C) const; + + // name sucks // Find the line which has as coordinate C pt[C] and gets it other coord // from planex. See if it is inside the sphere at t @@ -239,6 +256,8 @@ struct Arrangement_of_spheres_traits_3 { private: + + // HDS hds_; // ick, this is to handle location of points which are not already there // -1, -2 are bl, tr inf corners