// Copyright (c) 2014 INRIA Sophia-Antipolis (France), INRIA Lorraine LORIA. // All rights reserved. // // This file is part of CGAL (www.cgal.org). // You can redistribute it and/or modify it under the terms of the GNU // General Public License as published by the Free Software Foundation, // either version 3 of the License, or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // Author(s) : Fernando de Goes, Pierre Alliez, Ivo Vigan, Clément Jamin #ifndef CGAL_RECONSTRUCTION_TRIANGULATION_2_H #define CGAL_RECONSTRUCTION_TRIANGULATION_2_H // local #include #include #include #include #include // CGAL #include #include #include #include // boost #include #include #include #include #include #include #include // STL #include #include #include #include #include #define CGAL_EPS 1e-15 namespace CGAL { namespace OTR_2 { /// \internal /// The Reconstruction_triangulation_2 class /// provides the reconstruction simplex as well as the transportation plan. /// - Each vertex stores a normal vector. /// - A vertex a Sample which got assigned to it by the transportation plan, /// well as the corresponding relocated Point (of type Traits_::Point_2). /// - In order to solve a linear system over the triangulation, a vertex may be constrained /// or not (i.e. may contribute to the right or left member of the linear system), /// and has a unique index. /// The vertex class must derive from Reconstruction_vertex_base_3. /// /// @param Traits_ The traits class /// @param Tds Model of TriangulationDataStructure_2. /// The vertex class must derive from Reconstruction_vertex_base_2. /// The face class must derive from Reconstruction_face_base_2. template< class Traits_, class Tds_ = Triangulation_data_structure_2< Reconstruction_vertex_base_2, Reconstruction_face_base_2 > > class Reconstruction_triangulation_2 : public Delaunay_triangulation_2 { public: typedef Delaunay_triangulation_2 Base; typedef typename Traits_::FT FT; typedef typename Traits_::Point_2 Point; typedef typename Traits_::Vector_2 Vector; typedef typename Traits_::Line_2 Line; typedef typename Traits_::Segment_2 Segment; typedef typename Traits_::Triangle_2 Triangle; typedef typename Base::Vertex Vertex; typedef typename Base::Vertex_handle Vertex_handle; typedef typename Base::Vertex_iterator Vertex_iterator; typedef typename Base::Vertex_circulator Vertex_circulator; typedef typename Base::Finite_vertices_iterator Finite_vertices_iterator; typedef typename Base::Edge Edge; typedef typename Base::Edge_circulator Edge_circulator; typedef typename Base::Finite_edges_iterator Finite_edges_iterator; typedef typename Base::Face_handle Face_handle; typedef typename Base::Face_circulator Face_circulator; typedef typename Base::Finite_faces_iterator Finite_faces_iterator; typedef std::map > Vertex_handle_map; typedef std::map > Face_handle_map; typedef std::set > Vertex_handle_set; typedef std::set > Edge_set; typedef std::vector Edge_vector; typedef OTR_2::Cost Cost_; typedef OTR_2::Sample Sample_; typedef std::vector Sample_vector; typedef typename Sample_vector::const_iterator Sample_vector_const_iterator; typedef OTR_2::Sample_with_priority PSample; typedef std::priority_queue, OTR_2::greater_priority > SQueue; typedef Reconstruction_edge_2 Rec_edge_2; using Base::geom_traits; typedef boost::multi_index_container< Rec_edge_2, boost::multi_index::indexed_by< // sort by Rec_edge_2::operator< boost::multi_index::ordered_unique< boost::multi_index::identity< Rec_edge_2 > > , // sort by Rec_edge_2::priority() boost::multi_index::ordered_non_unique< boost::multi_index::const_mem_fun< Rec_edge_2,const FT,&Rec_edge_2::priority> > > > MultiIndex; FT m_factor; // ghost vs solid public: Reconstruction_triangulation_2(Traits_ traits = Traits_()) : Base(traits), m_factor(1.) { } ~Reconstruction_triangulation_2() { } FT& ghost_factor() { return m_factor; } FT ghost_factor() const { return m_factor; } Edge random_finite_edge() const { static CGAL::Random rng; std::size_t nbf = Base::number_of_faces(); int offset = rng.get_int(0, static_cast(nbf - 1)); Finite_faces_iterator fi = Base::finite_faces_begin(); for (int i = 0; i < offset; i++) fi++; Face_handle face = fi; int index = rng.get_int(0, 40) % 3; return Edge(face, index); } // ACCESS // Vertex_handle source_vertex(const Edge& edge) const { return edge.first->vertex(Base::ccw(edge.second)); } Vertex_handle target_vertex(const Edge& edge) const { return edge.first->vertex(Base::cw(edge.second)); } Vertex_handle opposite_vertex(const Edge& edge) const { return edge.first->vertex(edge.second); } bool is_pinned(const Edge& edge) const { Vertex_handle s = source_vertex(edge); if (s->pinned()) return true; return false; } Edge twin_edge(const Edge& edge) const { Face_handle f = edge.first; Vertex_handle v = source_vertex(edge); Face_handle nf = f->neighbor(edge.second); return Edge(nf, Base::ccw(nf->index(v))); } Edge next_edge(const Edge& edge) const { Face_handle f = edge.first; int index = Base::ccw(edge.second); return Edge(f, index); } Edge prev_edge(const Edge& edge) const { Face_handle f = edge.first; int index = Base::cw(edge.second); return Edge(f, index); } FT get_length(const Edge& edge) const { Segment segment = get_segment(edge); return CGAL::approximate_sqrt(segment.squared_length()); } Segment get_segment(const Edge& edge) const { const Point& ps = source_vertex(edge)->point(); const Point& pt = target_vertex(edge)->point(); return geom_traits().construct_segment_2_object()(ps, pt); } Triangle get_triangle(Face_handle face) const { Vertex_handle v0 = face->vertex(0); Vertex_handle v1 = face->vertex(1); Vertex_handle v2 = face->vertex(2); return geom_traits().construct_triangle_2_object()( v0->point(), v1->point(), v2->point()); } // GET LINK // void get_vertices_from_edge_link(const Edge& edge, Vertex_handle_set& vertices) const { vertices.insert(opposite_vertex(edge)); vertices.insert(opposite_vertex(twin_edge(edge))); } void get_vertices_from_vertex_link(Vertex_handle vertex, Vertex_handle_set& vertices) const { Vertex_circulator vcirc = Base::incident_vertices(vertex); Vertex_circulator vend = vcirc; CGAL_For_all(vcirc, vend) { Vertex_handle v = vcirc; vertices.insert(v); } } // boundary of star(vertex) // 'outward' chooses the orientation of the boundary void get_edges_from_star_minus_link(Vertex_handle vertex, Edge_vector& hull, bool outward = false) const { Face_circulator fcirc = Base::incident_faces(vertex); Face_circulator fend = fcirc; CGAL_For_all(fcirc, fend) { Face_handle face = fcirc; int index = face->index(vertex); Edge edge(face, index); if (outward) edge = twin_edge(edge); hull.push_back(edge); } } // ATTRIBUTES // bool is_ghost(const Edge& edge) const { return edge.first->ghost(edge.second); } int get_plan(const Edge& edge) const { return edge.first->plan(edge.second); } void set_plan(const Edge& edge, int simplex) const { edge.first->plan(edge.second) = simplex; } FT get_mass(const Edge& edge) const { return edge.first->mass(edge.second); } void set_mass(const Edge& edge, const FT mass) const { edge.first->mass(edge.second) = mass; } const Cost_& get_cost(const Edge& edge) const { return edge.first->cost(edge.second); } void set_vertex_cost(const Edge& edge, const Cost_& cost) const { edge.first->vertex_cost(edge.second) = cost; } void set_edge_cost(const Edge& edge, const Cost_& cost) const { edge.first->edge_cost(edge.second) = cost; } FT get_vertex_minus_edge_cost(const Edge& edge) const { const Cost_& vcost = edge.first->vertex_cost(edge.second); const Cost_& ecost = edge.first->edge_cost(edge.second); return vcost.finalize() - m_factor * ecost.finalize(); } FT get_vertex_over_edge_cost(const Edge& edge) const { FT vvalue = edge.first->vertex_cost(edge.second).finalize(); FT evalue = edge.first->edge_cost(edge.second).finalize(); if (evalue == vvalue) return FT(1) / m_factor; return vvalue / (m_factor * evalue); } FT get_edge_relevance(const Edge& edge) const { FT M = get_mass(edge); if (M == FT(0)) return FT(0); FT L = get_length(edge); FT cost = get_cost(edge).finalize(); return M * L * L / cost; } FT get_density(const Edge& edge) const { FT length = get_length(edge); FT mass = get_mass(edge); return (mass / length); } unsigned int nb_samples(const Edge& edge) const { Edge twin = twin_edge(edge); return edge.first->samples(edge.second).size() + twin.first->samples(twin.second).size(); } void collect_samples_from_edge( const Edge& edge, Sample_vector& samples) const { const Sample_vector& edge_samples = edge.first->samples(edge.second); samples.insert(samples.end(), edge_samples.begin(), edge_samples.end()); } void collect_samples_from_vertex( Vertex_handle vertex, Sample_vector& samples, bool cleanup) const { Face_circulator fcirc = Base::incident_faces(vertex); Face_circulator fend = fcirc; CGAL_For_all(fcirc, fend) { Face_handle face = fcirc; int index = face->index(vertex); Edge edge(face, index); collect_samples_from_edge(edge, samples); Edge next = next_edge(edge); collect_samples_from_edge(next, samples); Edge prev = prev_edge(edge); collect_samples_from_edge(prev, samples); if (cleanup) face->clean_all_samples(); } Sample_* sample = vertex->sample(); if (sample) samples.push_back(sample); if (cleanup) vertex->set_sample(NULL); } void collect_all_samples(Sample_vector& samples) const { for (Finite_edges_iterator ei = Base::finite_edges_begin(); ei != Base::finite_edges_end(); ++ei) { Edge edge = *ei; Edge twin = twin_edge(edge); collect_samples_from_edge(edge, samples); collect_samples_from_edge(twin, samples); } for (Finite_vertices_iterator vi = Base::finite_vertices_begin(); vi != Base::finite_vertices_end(); ++vi) { Vertex_handle v = vi; Sample_* sample = v->sample(); if (sample) samples.push_back(sample); } } void cleanup_assignments() { for (Finite_faces_iterator fi = Base::finite_faces_begin(); fi != Base::finite_faces_end(); ++fi) { fi->clean_all_samples(); } for (Finite_vertices_iterator vi = Base::finite_vertices_begin(); vi != Base::finite_vertices_end(); ++vi) { vi->set_sample(NULL); } } // COST // Cost_ compute_total_cost() const { Cost_ sum; for (Finite_edges_iterator ei = Base::finite_edges_begin(); ei != Base::finite_edges_end(); ++ei) { Edge edge = *ei; const Cost_& cost = get_cost(edge); sum.update_max(cost); sum.add(cost); } return sum; } Cost_ compute_cost_around_vertex(Vertex_handle vertex) const { Cost_ inner; Cost_ outer; Face_circulator fcirc = Base::incident_faces(vertex); Face_circulator fend = fcirc; CGAL_For_all(fcirc, fend) { Face_handle face = fcirc; int index = face->index(vertex); Edge edge(face, index); Cost_ cost = get_cost(edge); outer.update_max(cost); outer.add(cost); edge = next_edge(edge); cost = get_cost(edge); inner.update_max(cost); inner.add(cost); edge = next_edge(edge); cost = get_cost(edge); inner.update_max(cost); inner.add(cost); } inner.divide(2.0); Cost_ sum; sum.add(inner); sum.add(outer); sum.update_max(inner); sum.update_max(outer); return sum; } void reset_all_costs() { for (Finite_edges_iterator ei = Base::finite_edges_begin(); ei != Base::finite_edges_end(); ++ei) { Edge edge = *ei; update_cost(edge); } } void update_cost(const Edge& edge) { compute_mass(edge); compute_edge_cost(edge); compute_vertex_cost(edge); select_plan(edge); } void compute_mass(const Edge& edge) const { FT mass = 0; typename Sample_vector::const_iterator it; const Sample_vector& samples0 = edge.first->samples(edge.second); for (it = samples0.begin(); it != samples0.end(); ++it) { Sample_* sample = *it; mass += sample->mass(); } Edge twin = twin_edge(edge); const Sample_vector& samples1 = twin.first->samples(twin.second); for (it = samples1.begin(); it != samples1.end(); ++it) { Sample_* sample = *it; mass += sample->mass(); } set_mass(edge, mass); set_mass(twin, mass); } void select_plan(const Edge& edge) const { // transport plan: // 0 - to vertex // 1 - to edge int plan = 0; FT diff = get_vertex_minus_edge_cost(edge); if (diff >= 0.0) plan = 1; Edge twin = twin_edge(edge); set_plan(edge, plan); set_plan(twin, plan); } void compute_edge_cost(const Edge& edge) const { SQueue squeue; FT M = get_mass(edge); FT L = get_length(edge); sort_samples_from_edge(edge, squeue); Cost_ cost = compute_cost_from_squeue(squeue, M, L); Edge twin = twin_edge(edge); set_edge_cost(edge, cost); set_edge_cost(twin, cost); } void sort_samples_from_edge(const Edge& edge, SQueue& squeue) const { typename Sample_vector::const_iterator it; const Sample_vector& samples0 = edge.first->samples(edge.second); for (it = samples0.begin(); it != samples0.end(); ++it) { Sample_* sample = *it; squeue.push(PSample(sample, sample->coordinate())); } Edge twin = twin_edge(edge); const Sample_vector& samples1 = twin.first->samples(twin.second); for (it = samples1.begin(); it != samples1.end(); ++it) { Sample_* sample = *it; squeue.push(PSample(sample, 1.0 - sample->coordinate())); } } Cost_ compute_cost_from_squeue(SQueue& squeue, const FT M, const FT L) const { if (squeue.empty()) return Cost_(); if (M == FT(0)) return Cost_(); Cost_ sum; FT start = 0; FT coef = L / M; while (!squeue.empty()) { PSample psample = squeue.top(); squeue.pop(); FT mass = psample.sample()->mass(); FT coord = psample.priority() * L; FT bin = mass * coef; FT center = start + FT(0.5) * bin; FT pos = coord - center; FT norm2 = psample.sample()->distance2(); FT tang2 = bin * bin / 12 + pos * pos; sum.add(Cost_(norm2, tang2), mass); sum.compute_max(norm2, tang2); start += bin; } return sum; } void compute_vertex_cost(const Edge& edge) const { Edge twin = twin_edge(edge); const Point& ps = source_vertex(edge)->point(); const Point& pt = target_vertex(edge)->point(); Sample_vector samples; collect_samples_from_edge(edge, samples); collect_samples_from_edge(twin, samples); Cost_ sum; for (Sample_vector_const_iterator it = samples.begin(); it != samples.end(); ++it) { Sample_* sample = *it; FT mass = sample->mass(); const Point& query = sample->point(); FT Ds = geom_traits().compute_squared_distance_2_object()(query, ps); FT Dt = geom_traits().compute_squared_distance_2_object()(query, pt); FT dist2 = ((std::min))(Ds, Dt); FT norm2 = sample->distance2(); FT tang2 = dist2 - norm2; sum.add(Cost_(norm2, tang2), mass); sum.compute_max(norm2, tang2); } set_vertex_cost(edge, sum); set_vertex_cost(twin, sum); } // SAMPLE // template // value_type = Sample_* void assign_samples(Iterator begin, Iterator end) { for (Iterator it = begin; it != end; ++it) { Sample_* sample = *it; assign_sample(sample); } } template // value_type = Sample_* void assign_samples_brute_force(Iterator begin, Iterator end) { for (Iterator it = begin; it != end; ++it) { Sample_* sample = *it; assign_sample_brute_force(sample); } } bool assign_sample(Sample_* sample) { const Point& point = sample->point(); Face_handle face = Base::locate(point); if (face == Face_handle() || Base::is_infinite(face)) { //std::cout << "free bird" << std::endl; return false; } Vertex_handle vertex = find_nearest_vertex(point, face); if (vertex != Vertex_handle()) { assign_sample_to_vertex(sample, vertex); return true; } Edge edge = find_nearest_edge(point, face); assign_sample_to_edge(sample, edge); return true; } bool assign_sample_brute_force(Sample_* sample) { const Point& point = sample->point(); Face_handle nearest_face = Face_handle(); for (Finite_faces_iterator fi = Base::finite_faces_begin(); fi != Base::finite_faces_end(); ++fi) { Face_handle face = fi; if (face_has_point(face, point)) { nearest_face = face; break; } } if (nearest_face == Face_handle()) { //std::cout << "free bird" << std::endl; return false; } Vertex_handle vertex = find_nearest_vertex(point, nearest_face); if (vertex != Vertex_handle()) { assign_sample_to_vertex(sample, vertex); return true; } Edge edge = find_nearest_edge(point, nearest_face); assign_sample_to_edge(sample, edge); return true; } bool face_has_point(Face_handle face, const Point& query) const { for (int i = 0; i < 3; ++i) { Edge edge(face, i); const Point& ps = source_vertex(edge)->point(); const Point& pt = target_vertex(edge)->point(); if (!compute_triangle_ccw(ps, pt, query)) return false; } return true; } Vertex_handle find_nearest_vertex(const Point& point, Face_handle face) const { for (int i = 0; i < 3; ++i) { Vertex_handle vi = face->vertex(i); const Point& pi = vi->point(); if (pi == point) return vi; } return Vertex_handle(); } Edge find_nearest_edge(const Point& point, Face_handle face) const { Edge nearest(face, 0); FT min_dist2 = compute_distance2(point, get_segment(nearest)); for (int i = 1; i < 3; ++i) { Edge edge(face, i); Segment segment = get_segment(edge); FT dist2 = compute_distance2(point, segment); if (dist2 < min_dist2) { min_dist2 = dist2; nearest = edge; } } return nearest; } void assign_sample_to_vertex(Sample_* sample, Vertex_handle vertex) const { /*if (vertex->sample()) { std::cout << "assign to vertex: vertex already has sample" << std::endl; }*/ sample->distance2() = FT(0); sample->coordinate() = FT(0); vertex->set_sample(sample); } void assign_sample_to_edge(Sample_* sample, const Edge& edge) const { Segment segment = get_segment(edge); const Point& query = sample->point(); sample->distance2() = compute_distance2(query, segment); sample->coordinate() = compute_coordinate(query, segment); edge.first->add_sample(edge.second, sample); } FT compute_distance2(const Point& query, const Segment& segment) const { if (geom_traits().orientation_2_object()(segment.source(), segment.target(), query) == COLLINEAR) return FT(0); Line line = geom_traits().construct_line_2_object()(segment); return geom_traits().compute_squared_distance_2_object()(query, line); } FT compute_coordinate(const Point& q, const Segment& segment) const { const Point& p0 = segment.source(); const Point& p1 = segment.target(); Vector p0p1 = geom_traits().construct_vector_2_object()(p0, p1); Vector p0q = geom_traits().construct_vector_2_object()(p0, q); FT t = geom_traits().compute_scalar_product_2_object()(p0q, p0p1) / geom_traits().compute_scalar_product_2_object()(p0p1, p0p1); return t; // [0,1] } // SIGNED DISTANCE // // signed distance from line(a,b) to point t FT signed_distance(Vertex_handle a, Vertex_handle b, Vertex_handle t) const { const Point& pa = a->point(); const Point& pb = b->point(); const Point& pt = t->point(); return compute_signed_distance(pa, pb, pt); } // signed distance from line(a,b) to point t FT compute_signed_distance( const Point& pa, const Point& pb, const Point& pt) const { if (pt == pa) return FT(0); if (pt == pb) return FT(0); if (pa == pb) return CGAL::approximate_sqrt(geom_traits().compute_squared_distance_2_object()(pa, pt)); Vector vab = geom_traits().construct_vector_2_object()(pa, pb); // Normalize vab vab = geom_traits().construct_scaled_vector_2_object()( vab, FT(1) / CGAL::approximate_sqrt(geom_traits().compute_squared_length_2_object()(vab))); Vector vab90 = geom_traits().construct_vector_2_object()(-vab.y(), vab.x()); Vector vat = geom_traits().construct_vector_2_object()(pa, pt); return geom_traits().compute_scalar_product_2_object()(vat, vab90); } // signed distance from t to the intersection of line(a,b) and line(t,s) // the pair::first is false if sign==-1 and true otherwise std::pair > signed_distance_from_intersection(Vertex_handle a, Vertex_handle b, Vertex_handle t, Vertex_handle s) const { const Point& pa = a->point(); const Point& pb = b->point(); const Point& pt = t->point(); const Point& ps = s->point(); return compute_signed_distance_from_intersection(pa, pb, pt, ps); } // signed distance from t to the intersection of line(a,b) and line(t,s) // the pair::first is false if sign==-1 and true otherwise std::pair > compute_signed_distance_from_intersection( const Point& pa, const Point& pb, const Point& pt, const Point& ps) const { FT Dabt = compute_signed_distance(pa, pb, pt); if (Dabt == FT(0)) return std::make_pair(true,boost::make_optional(FT(0))); Line lab = geom_traits().construct_line_2_object()( pa, geom_traits().construct_vector_2_object()(pa, pb)); Line lts = geom_traits().construct_line_2_object()( pt, geom_traits().construct_vector_2_object()(pt, ps)); boost::optional Dqt; typename CGAL::cpp11::result_of::type result = intersection(lab, lts); if (result) { const Point* iq = boost::get(&(*result)); if (iq) Dqt = CGAL::approximate_sqrt(geom_traits().compute_squared_distance_2_object()(*iq, pt)); } return std::make_pair( (Dabt < FT(0) ? false : true) ,Dqt ); } bool is_triangle_ccw(Vertex_handle a, Vertex_handle b, Vertex_handle c) const { const Point& pa = a->point(); const Point& pb = b->point(); const Point& pc = c->point(); return compute_triangle_ccw(pa, pb, pc); } bool compute_triangle_ccw( const Point& pa, const Point& pb, const Point& pc) const { return geom_traits().orientation_2_object()(pa,pb,pc) != RIGHT_TURN; } // COMBINATORIAL TESTS // // (a,b) is cyclic if (a,b,c) and (a,c,b) exist bool is_edge_cyclic(const Edge& edge) const { Vertex_handle f = opposite_vertex(edge); Vertex_handle b = opposite_vertex(twin_edge(edge)); return (f == b); } // b from (a,b) is cyclic if (a,b,c) and (b,a,c) exist bool is_target_cyclic(const Edge& edge) const { if (!is_edge_cyclic(edge)) return false; Edge twin = twin_edge(edge); Edge prev = prev_edge(twin); Face_handle fp = prev.first->neighbor(prev.second); Face_handle ft = twin.first->neighbor(twin.second); return (fp == ft); } bool is_flippable(const Edge& edge) const { Edge twin = twin_edge(edge); if (Base::is_infinite(twin.first)) return false; if (Base::is_infinite(edge.first)) return false; Vertex_handle vs = source_vertex(edge); Vertex_handle vt = target_vertex(edge); Vertex_handle vf = opposite_vertex(edge); Vertex_handle vb = opposite_vertex(twin); return is_triangle_ccw(vs, vb, vf) && is_triangle_ccw(vt, vf, vb); } bool is_collapsible(const Edge& edge) const { return check_link_test(edge) && check_kernel_test(edge); } bool check_link_test(const Edge& edge) const { Vertex_handle s = source_vertex(edge); Vertex_handle t = target_vertex(edge); if (s == t) return false; typename Vertex_handle_set::const_iterator it; Vertex_handle_set svertices; get_vertices_from_vertex_link(s, svertices); Vertex_handle_set tvertices; get_vertices_from_vertex_link(t, tvertices); // link(s) inter link(t) Vertex_handle_set ivertices; for (it = svertices.begin(); it != svertices.end(); ++it) { Vertex_handle v = *it; if (tvertices.find(v) != tvertices.end()) ivertices.insert(v); } Vertex_handle_set evertices; get_vertices_from_edge_link(edge, evertices); // link(edge) =? link(s) inter link(t) if (evertices.size() != ivertices.size()) return false; for (it = evertices.begin(); it != evertices.end(); ++it) { Vertex_handle v = *it; if (ivertices.find(v) == ivertices.end()) return false; } return true; } bool check_kernel_test(const Edge& edge) const { Vertex_handle s = source_vertex(edge); Vertex_handle t = target_vertex(edge); Edge_vector hull; hull.reserve(16); get_edges_from_star_minus_link(s, hull); return is_in_kernel(t->point(), hull.begin(), hull.end()); } template // value_type = Edge bool is_in_kernel(const Point& query, Iterator begin, Iterator end) const { for (Iterator it = begin; it != end; ++it) { Edge edge = *it; const Point& pa = source_vertex(edge)->point(); const Point& pb = target_vertex(edge)->point(); if (!compute_triangle_ccw(pa, pb, query)) return false; } return true; } // COLLAPSE // // (s,a,b) + (s,b,c) -> (s,a,c) + (a,b,c) // st = (source,target) from 'make_collapsible' // return (a,c) Edge flip(const Edge& sb, Edge& st, int /*verbose*/ = 0) { Vertex_handle t = target_vertex(st); Edge sc = twin_edge(prev_edge(sb)); Base::tds().flip(sb.first, sb.second); Edge ac = prev_edge(twin_edge(sc)); Vertex_handle a = source_vertex(ac); if (a == t) st = prev_edge(ac); return ac; } void collapse(const Edge& edge, int /*verbose*/ = 0) { if (is_edge_cyclic(edge)) { collapse_cyclic_edge(edge); return; } Edge twin = twin_edge(edge); Base::tds().join_vertices(twin); } // (a,b,c) + (c,b,a) + (a,c,i) + (c,a,j) -> // (a,c,i) + (c,a,j) void collapse_cyclic_edge(const Edge& bc, int verbose = 0) { if (verbose > 1) std::cout << "collapse_cyclic_edge ... "; Edge cb = twin_edge(bc); Face_handle abc = bc.first; Face_handle cba = cb.first; Vertex_handle b = source_vertex(bc); Vertex_handle c = target_vertex(bc); Vertex_handle a = opposite_vertex(bc); Edge ac = twin_edge(next_edge(bc)); Edge ca = twin_edge(prev_edge(cb)); a->set_face(ac.first); c->set_face(ca.first); ac.first->set_neighbor(ac.second, ca.first); ca.first->set_neighbor(ca.second, ac.first); this->delete_face(abc); this->delete_face(cba); this->delete_vertex(b); if (verbose > 1) std::cout << "done" << std::endl; } void print_edge(Rec_edge_2 edge) const { int i = ((edge).edge()).second; Point a = ((edge).edge()).first->vertex((i+1)%3)->point(); Point b = ((edge).edge()).first->vertex((i+2)%3)->point(); std::cout <<"( " << (edge).priority() << ") ( " << a << " , " << b << " )" << std::endl; } bool is_p_infinity(const std::pair >& p) const { return p.first && p.second==boost::none; } bool is_m_infinity(const std::pair >& p) const { return !p.first && p.second==boost::none; } bool is_infinity(const std::pair >& p) const { return p.second==boost::none; } std::pair > m_infinity() const { return std::pair >(false,boost::optional()); } template // value_type = Edge bool make_collapsible(Edge& edge, Iterator begin, Iterator end, int verbose = 0) { Vertex_handle source = source_vertex(edge); Vertex_handle target = target_vertex(edge); MultiIndex multi_ind; for (Iterator it = begin; it != end; ++it) { Edge ab = twin_edge(*it); Vertex_handle a = source_vertex(ab); Vertex_handle b = target_vertex(ab); std::pair > D = signed_distance_from_intersection(a, b, target, source); if (!D.first ) { CGAL_assertion(D.second!=boost::none); multi_ind.insert(Rec_edge_2(ab, *D.second)); } } int nb_flips = 0; while (!multi_ind.empty()) { Rec_edge_2 pedge = *(multi_ind.template get<1>()).begin(); FT Dbc = pedge.priority(); Edge bc = pedge.edge(); (multi_ind.template get<0>()).erase(pedge); Edge sb = prev_edge(bc); Edge ab = prev_edge(twin_edge(sb)); Edge sc = twin_edge(next_edge(bc)); Edge cd = next_edge(sc); Vertex_handle a = source_vertex(ab); Vertex_handle b = source_vertex(bc); Vertex_handle c = target_vertex(bc); Vertex_handle d = target_vertex(cd); std::pair > Dac=m_infinity(); if (a != c && is_triangle_ccw(a, b, c)) Dac = signed_distance_from_intersection(a, c, target, source); std::pair > Dbd=m_infinity(); if (b != d && is_triangle_ccw(b, c, d)) Dbd = signed_distance_from_intersection(b, d, target, source); if ( is_m_infinity(Dac) && is_m_infinity(Dbd) ) { if (verbose > 1) std::cerr << "--- No flips available ---" << std::endl; return false; } FT value = Dbc <= 0 ? 1 : 2*Dbc; // value used if Dbd or Dac are +infinity if ( !is_infinity(Dac) ) { if ( !is_infinity(Dbd)) value = (std::max)(*Dac.second * (Dac.first?1:-1), *Dbd.second * (Dbd.first?1:-1) ); else if ( !is_p_infinity(Dbd) ) value = *Dac.second * (Dac.first?1:-1); } else if ( !is_infinity(Dbd) && !is_p_infinity(Dac)) value = *Dbd.second * (Dbd.first?1:-1); // if ( max(Dac,Dbd)+CGAL_EPS < Dbc ) if (value + CGAL_EPS < Dbc) { /* std::cerr.precision(10); std::cerr << "--- Flip makes kernel worse ---" << std::endl; std::cerr << Dac << " or " << Dbd << " vs " << Dbc << std::endl; std::cerr << "a: " << a->point() << std::endl; std::cerr << "b: " << b->point() << std::endl; std::cerr << "c: " << c->point() << std::endl; std::cerr << "d: " << d->point() << std::endl; std::cerr << "t: " << target->point() << std::endl; std::cerr << "diff = " << Dbc - (std::max)(Dac, Dbd) << std::endl; */ return false; } // if (Dac > Dbd) if (is_p_infinity(Dac) || is_m_infinity(Dbd) || (!is_p_infinity(Dbd) && !is_m_infinity(Dac) && *Dac.second * (Dac.first?1:-1) > *Dbd.second * (Dbd.first?1:-1))) { (multi_ind.template get<0>()).erase(Rec_edge_2(ab)); Edge ac = flip(sb, edge, verbose); if (!Dac.first) { multi_ind.insert(Rec_edge_2(ac, - *Dac.second)); } } else { (multi_ind.template get<0>()).erase(Rec_edge_2(cd)); Edge bd = flip(sc, edge, verbose); if (!Dbd.first) { multi_ind.insert(Rec_edge_2(bd, - *Dbd.second)); } } nb_flips++; } if (verbose > 1) std::cerr << "Nb flips: " << nb_flips << std::endl; return true; } }; } } //namespace CGAL #undef CGAL_EPS #endif // CGAL_RECONSTRUCTION_TRIANGULATION_2_H