From 68daecef6e6b23a17a5d159399ec5446bc9d4017 Mon Sep 17 00:00:00 2001 From: Thijs van Lankveld Date: Mon, 25 Feb 2013 18:26:51 +0100 Subject: [PATCH] Very first version of conforming triangulation 3. --- .../Conforming_Delaunay_triangulation_3.h | 1264 +++++++++++++++++ .../Conforming_triangulation_cell_base_3.h | 179 +++ .../Conforming_triangulation_vertex_base_3.h | 69 + .../CGAL/Delaunay_triangulation_utils_3.h | 402 ++++++ .../include/CGAL/Encroaching_collecter_3.h | 123 ++ .../CGAL/Triangulation_segment_traverser_3.h | 695 +++++++++ .../Conforming_triangulation_3/copyright | 2 + .../description.txt | 3 + .../Conforming_triangulation_3/license.txt | 1 + .../Conforming_triangulation_3/maintainer | 1 + 10 files changed, 2739 insertions(+) create mode 100644 Conforming_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h create mode 100644 Conforming_triangulation_3/include/CGAL/Conforming_triangulation_cell_base_3.h create mode 100644 Conforming_triangulation_3/include/CGAL/Conforming_triangulation_vertex_base_3.h create mode 100644 Conforming_triangulation_3/include/CGAL/Delaunay_triangulation_utils_3.h create mode 100644 Conforming_triangulation_3/include/CGAL/Encroaching_collecter_3.h create mode 100644 Conforming_triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h create mode 100644 Conforming_triangulation_3/package_info/Conforming_triangulation_3/copyright create mode 100644 Conforming_triangulation_3/package_info/Conforming_triangulation_3/description.txt create mode 100644 Conforming_triangulation_3/package_info/Conforming_triangulation_3/license.txt create mode 100644 Conforming_triangulation_3/package_info/Conforming_triangulation_3/maintainer diff --git a/Conforming_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h b/Conforming_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h new file mode 100644 index 00000000000..c981418b4fe --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h @@ -0,0 +1,1264 @@ +//The conforming Delaunay triangulation structure. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +// Based on CGAL/Delaunay_triangulation_3.h +// +// For the love of something, why is the 3D triangulation package so badly written? + +#ifndef CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H +#define CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H + +#include + +#include +#include + +#include "Delaunay_triangulation_utils_3.h" +#include "Conforming_triangulation_cell_base_3.h" +#include "Conforming_triangulation_vertex_base_3.h" +#include "Encroaching_collecter_3.h" + +CGAL_BEGIN_NAMESPACE + +#ifndef CGAL_CONSTRAINED_TRIANGULATION_2_H +struct No_intersection_tag{}; +struct Exact_intersections_tag{}; // To be used with an exact number type. +struct Exact_predicates_tag{}; // To be used with filtered exact number type. +#endif + +template < class Tr > class Natural_neighbors_3; + +template < class Gt, + class Tds = Triangulation_data_structure_3 < Conforming_triangulation_vertex_base_3, + Conforming_triangulation_cell_base_3 >, + class Itag = No_intersection_tag > +class Conforming_Delaunay_triangulation_3: public Delaunay_triangulation_utils_3 { + typedef Triangulation_data_structure Tds; + + typedef Conforming_Delaunay_triangulation_3 cDT; + typedef Delaunay_triangulation_utils_3 DT; + typedef Triangulation_3 Tr; + + friend class Natural_neighbors_3; + +public: + typedef Tds Triangulation_data_structure; + typedef Gt Geom_traits; + typedef Itag Intersection_tag; + + typedef typename Gt::FT FT; + + typedef typename Gt::Point_3 Point; + typedef typename Gt::Vector_3 Vector; + typedef typename Gt::Segment_3 Segment; + typedef typename Gt::Triangle_3 Triangle; + typedef typename Gt::Tetrahedron_3 Tetrahedron; + typedef typename Gt::Plane_3 Plane; + + typedef typename DT::Cell_handle Cell_handle; + typedef typename DT::Vertex_handle Vertex_handle; + + typedef typename DT::Cell Cell; + typedef typename DT::Vertex Vertex; + typedef typename DT::Facet Facet; + typedef typename DT::Edge Edge; + + typedef typename DT::Bi_vertex Bi_vertex; + + typedef std::pair Constraint_1; + + typedef typename DT::Cell_circulator Cell_circulator; + typedef typename DT::Facet_circulator Facet_circulator; + typedef typename DT::Cell_iterator Cell_iterator; + typedef typename DT::Facet_iterator Facet_iterator; + typedef typename DT::Edge_iterator Edge_iterator; + typedef typename DT::Vertex_iterator Vertex_iterator; + + typedef typename DT::Finite_vertices_iterator Finite_vertices_iterator; + typedef typename DT::Finite_cells_iterator Finite_cells_iterator; + typedef typename DT::Finite_facets_iterator Finite_facets_iterator; + typedef typename DT::Finite_edges_iterator Finite_edges_iterator; + + typedef typename DT::All_cells_iterator All_cells_iterator; + + typedef typename DT::Cell_traverser Cell_traverser; + + typedef typename DT::Locate_type Locate_type; + +private: + typedef std::list Point_list; + typedef typename Point_list::iterator Points_iterator; + typedef std::back_insert_iterator Points_output; + +public: + typedef Encroaching_collecter_3 Encroaching_collecter; + +#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 + using DT::cw; + using DT::ccw; + using DT::geom_traits; + using DT::number_of_vertices; + using DT::dimension; + using DT::finite_facets_begin; + using DT::finite_facets_end; + using DT::finite_vertices_begin; + using DT::finite_vertices_end; + using DT::finite_cells_begin; + using DT::finite_cells_end; + using DT::finite_edges_begin; + using DT::finite_edges_end; + using DT::tds; + using DT::infinite_vertex; + using DT::insert; + using DT::next_around_edge; + using DT::vertex_triple_index; + using DT::mirror_vertex; + using DT::move_point; + using DT::coplanar; + using DT::coplanar_orientation; + using DT::orientation; +#endif + +protected: + // Visitor that should be run before and after inserting a point + // to maintain the conforming edges in the triangulation. + template < class Tr > + class Conforming_visitor { + typedef std::pair CBV; + + mutable std::list _conforming; + + protected: + Tr *_tr; + + public: + Conforming_visitor(Tr* tr): _tr(tr) {} + + // Store and remove the conforming edges from a collection of cells. + template + void process_cells_in_conflict(InputIterator start, InputIterator end) const { + for (InputIterator it = start; it != end; ++it) { + for (int i = 0; i < _tr->dimension(); ++i) { + for (int j = i+1; j < _tr->dimension()+1; ++j) { + if ((*it)->is_conforming(i, j)) { + _conforming.push_back(CBV(_tr->to_bi_vertex(*it, i, j), (*it)->is_marked(i, j))); + _tr->set_conforming(*it, i, j, false); + } + } + } + } + } + + // Reinsert the conforming edges into the triangulation. + void reinsert_vertices(Vertex_handle) { + while (!_conforming.empty()) { + CBV b = _conforming.front(); + _conforming.pop_front(); + if (b.second) + _tr->insert_marked(b.first); + else + _tr->insert_conforming(b.first); + } + } + + // Reinsert the conforming edges into the triangulation + // under the assumption that no Steiner points need to be added. + void reinsert_no_steiner() { + Cell_handle c; int li, lj; + while (!_conforming.empty()) { + CBV b = _conforming.front(); + _conforming.pop_front(); + CGAL_triangulation_assertion_code(bool ok =) + _tr->is_edge(b.first.first, b.first.second, c, li, lj); + CGAL_triangulation_assertion(ok); + + if (b.second) + _tr->mark_edge(c, li, lj); + else + _tr->set_conforming(c, li, lj, true); + } + } + + Vertex_handle replace_vertex(Cell_handle c, int index, const Point&) {return c->vertex(index);} + void hide_point(Cell_handle, const Point&) {} + }; // class Conforming_visitor + + friend class Conflict_tester_3; + friend class Conflict_tester_2; + friend class Encroaching_collecter; + + Conforming_visitor conforming_visitor; + +protected: + // PREDICATES + + // Checks if point c is encroaching upon the segment ab. + // A point is encroaching if it is inside the ball with diameter ab. + bool is_encroaching(const Point& a, const Point& b, const Point& c) const { + CGAL_triangulation_assertion(a != b); + CGAL_triangulation_assertion(a != c); + CGAL_triangulation_assertion(b != c); + return side_of_bounded_sphere(a, b, c) != ON_UNBOUNDED_SIDE; + } + +public: + // CONSTRUCTORS + Conforming_Delaunay_triangulation_3(const Gt& gt = Gt()): DT(gt), conforming_visitor(this) {} + + // Copy constructor duplicates vertices and cells. + Conforming_Delaunay_triangulation_3(const Conforming_Delaunay_triangulation_3& tr): DT(tr), conforming_visitor(this) { + CGAL_triangulation_postcondition(is_valid()); + } + + // Create a conforming 3D Delaunay triangulation from a number of constraints. + template < class InputIterator > + Conforming_Delaunay_triangulation_3(InputIterator begin, InputIterator end, const Gt& gt = Gt()): + DT(gt), conforming_visitor(this) { + insert_conforming(first, last); + CGAL_triangulation_postcondition(is_valid()); + } + +protected: + // Compute the intersection of ab with a line or plane. + bool compute_intersection(const Line& ab, const Line& l, Point& p) const { + Object result = Gt().intersect_3_object()(ab, l); + return assign(p, result); + } + bool compute_intersection(const Line& ab, const Plane& pl, Point& p) const { + Object result = Gt().intersect_3_object()(ab, pl); + return assign(p, result); + } + + // Compute the intersection of ab with either a line (edge) or plane (facet). + bool compute_intersection(Locate_type lt, Cell_handle c, int li, int lj, Vertex_handle va, Vertex_handle vb, Point& p) const { + if (lt == EDGE) + return compute_intersection(Line(va->point(), vb->point()), + Line(c->vertex(li)->point(), c->vertex(lj)->point()), p); + else + return compute_intersection(Line(va->point(), vb->point()), plane(c, li), p); + } + + // Insert the intersection of segment ab with either a facet or an edge. + Vertex_handle intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb) { + return intersect(c, lt, li, lj, va, vb, Itag());} + Vertex_handle intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, No_intersection_tag); + Vertex_handle intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, Exact_intersections_tag); + Vertex_handle intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, Exact_predicates_tag); + + // Split segment ab based on a reference point. + Point construct_split_point(Vertex_handle va, Vertex_handle vb, const Point& ref); + +protected: + // Conform an edge to important encroaching points. + void conform_segment(Vertex_handle va, Vertex_handle vb, Points_iterator encr_begin, Points_iterator encr_end); + + // Mark a segment that it is assumed to already be in the triangulation. + void mark_segment(Vertex_handle va, Vertex_handle vb); + // Insert a marked segment; + // this includes conforming the segment and marking its subsegments as conforming. + void insert_marked(Vertex_handle va, Vertex_handle vb); + void insert_marked(const Bi_vertex& m) {insert_marked(m.first, m.second);} + + // Mark an existing edge or facet as (un)constrained. + void set_conforming(Cell_handle c, int li, int lj, bool C); + void set_conforming(const Edge& e, bool C) {set_conforming(e.first, e.second, e.third, C);} + void mark_edge(Cell_handle c, int li, int lj); + void mark_edge(const Edge& e) {mark_edge(e.first, e.second, e.third);} + +protected: + // Restore the conforming edges of a facet based on the neighboring cell. + void restore_conforming(Cell_handle c, int li); + void restore_conforming(const Facet& f) {restore_conforming(f.first, f.second);} + +public: + // INSERT POINT + + virtual Vertex_handle insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj); + +public: + // INSERT CONFORMING + + // Insert a conforming segment. + void insert_conforming(Vertex_handle va, Vertex_handle vb) { + Point_list encr; conform_segment(va, vb, encr.begin(), encr.end());} + void insert_conforming(const Bi_vertex& c) { + insert_conforming(c.first, c.second);} + void insert_conforming(const Point& a, const Point& b, Cell_handle hint = Cell_handle()); + void insert_conforming(const Constraint_1& c, Cell_handle hint = Cell_handle()) { + insert_conforming(c.first, c.second, hint);} + + // Insert a number of conforming segments, presented as a collection of pairs of Point_3. + template < class InputIterator > + int insert_conforming(InputIterator begin, InputIterator end, Cell_handle hint = Cell_handle()); + + // Insert a loop of conforming segments, presented as a collection of Point_3. + template < class InputIterator > + int insert_conforming_loop(InputIterator begin, InputIterator end, Cell_handle hint = Cell_handle()); + template < class InputIterator > + int insert_marked_loop(InputIterator begin, InputIterator end, Cell_handle hint = Cell_handle()); + template < class InputIterator > + int remove_conforming_loop(InputIterator begin, InputIterator end, Cell_handle hint = Cell_handle()); + +protected: + template < class InputIterator, class OutputIterator > + void insert_loop_points(InputIterator begin, InputIterator end, OutputIterator out, Cell_handle hint = Cell_handle()); + +public: + // REMOVE + + // Vertices that are incident to a conforming edge cannot be removed. + bool remove(Vertex_handle v); + + // Remove a conforming edge or segment. + bool remove_conforming(Cell_handle c, int li, int lj); + bool remove_conforming(const Edge& e) { + return remove_conforming(e.first, e.second, e.third);} + bool remove_conforming(Vertex_handle va, Vertex_handle vb); + bool remove_conforming(Bi_vertex c) { + return remove_conforming(c.first, c.second);} + +public: + // QUERY CONFORMING + + // Check if an edge is conforming. + bool is_conforming(Cell_handle c, int li, int lj) const {return c->is_conforming(li, lj);} + bool is_conforming(const Edge& e) const {return is_conforming(e.first, e.second, e.third);} + + // Check if an edge is marked. + bool is_marked(Cell_handle c, int li, int lj) const {return c->is_marked(li, lj);} + bool is_marked(const Edge& e) const {return is_marked(e.first, e.second, e.third);} + + // Check if a vertex is incident to a conforming edge. + bool are_there_incident_conforming(Vertex_handle v) const; + + // Check if the triangulation contains an edge ac on the segment ab, where c(li, lj) is that edge. + bool includes_edge(Vertex_handle va, Vertex_handle vb, Vertex_handle& vc, Cell_handle& c, int& li, int& lj) const; + + // Check if the triangulation is valid. + virtual bool is_valid(bool verbose = false, int level = 0) const; + virtual bool is_valid(Cell_handle c, bool verbose = false, int level = 0) const; + +private: + // Collect the points encroaching upon the segment ac, where c is the first point on the ray ab. + // encroaching must point to a Point. + // Returns true if a conforming edge or constrained facet was encountered and the intersection inserted. + template + bool collect_encroaching(Vertex_handle va, Vertex_handle vb, OutputIterator encroaching, Vertex_handle& vc); + + // Sort the points encroaching upon the segment ab to the from of the collection. + // begin and end must point to Point. + Points_iterator sort_encroaching(Points_iterator begin, Points_iterator end, const Point& a, const Point& b) const; +}; // Conforming_Delaunay_triangulation_3 + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Vertex_handle +Conforming_Delaunay_triangulation_3:: +intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, No_intersection_tag) { + std::cerr << " sorry, this triangulation does not deal with" << std::endl + << " intersecting conforming edges" << std::endl; + CGAL_triangulation_assertion(false); + return Vertex_handle(); +} + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Vertex_handle +Conforming_Delaunay_triangulation_3:: +intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, Exact_intersections_tag) { + // Compute the intersection of the segment ab and a facet or edge of cell c. + // Note that no further check is made whether this intersection actually occurs inside the segment(s) or triangle. + Point p; + CGAL_triangulation_assertion_code( bool ok = ) + compute_intersection(lt, c, li, lj, va, vb, p); + CGAL_triangulation_assertion(ok); + Vertex_handle v = insert(p, c); + v->steiner() = true; + return v; +} + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Vertex_handle +Conforming_Delaunay_triangulation_3:: +intersect(Cell_handle c, Locate_type lt, int li, int lj, Vertex_handle va, Vertex_handle vb, Exact_predicates_tag) { + // Compute the intersection of the segment ab and a facet or edge of cell c. + // Note that no further check is made whether this intersection actually occurs inside the segment(s) or triangle. + Point p; + bool ok = compute_intersection(lt, c, li, lj, va, vb, p); + if (ok) { + Vertex_handle v = insert(p, c); + v->steiner() = true; + return v; + } + + // Return the vertex closest to the intersection. + Vertex_handle closest; + switch (lt) { + case EDGE: { + Line ab(va->point(), vb->point()); + Line ij(c->vertex(li)->point(), c->vertex(lj)->point()); + FT dist2 = squared_distance(va->point(), ij); + FT db = squared_distance(vb->point(), ij); + FT di = squared_distance(c->vertex(li)->point(), ab); + FT dj = squared_distance(c->vertex(lj)->point(), ab); + closest = va; + if (db < dist2) {dist2 = db; closest = vb;} + if (di < dist2) {dist2 = di; closest = c->vertex(li);} + if (dj < dist2) {dist2 = dj; closest = c->vertex(lj);} + } + case FACET: { + Line ab(va->point(), vb->point()); + Plane p = plane(c, li); + FT dist2 = squared_distance(va->point(), p); + FT db = squared_distance(vb->point(), p); + FT d1 = squared_distance(c->vertex((li+1)&3)->point(), ab); + FT d2 = squared_distance(c->vertex((li+2)&3)->point(), ab); + FT d3 = squared_distance(c->vertex((li+3)&3)->point(), ab); + closest = va; + if (db < dist2) {dist2 = db; closest = vb;} + if (d1 < dist2) {dist2 = d1; closest = c->vertex((li+1)&3);} + if (d2 < dist2) {dist2 = d2; closest = c->vertex((li+2)&3);} + if (d3 < dist2) {dist2 = d3; closest = c->vertex((li+3)&3);} + } + default: + CGAL_triangulation_assertion(false); + } + + closest->steiner() = true; + return closest; +} + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Point +Conforming_Delaunay_triangulation_3:: +construct_split_point(Vertex_handle va, Vertex_handle vb, const Point& ref) { + // Blend of Si's and Shewchuk's segment recovery. + // We are actually interested in constructing a layer of protecting Steiner points around each acute vertex, + // such that no new Steiner point is inserted inside this layer. + // Shewchuk and Si place the protecting Steiner points on a ball around the acute vertex. + // However, determining if a vertex is acute is expensive. Additionally, the main purpose of the protecting points + // is to place new points outside the diametral ball with their acute vertex. + // This is always done if we place new Steiner points on the projection of the reference point onto the segment. + // To guarantee a minimum edge length, we never place a Steiner point closer to an endpoint than to the reference. + +#ifdef SPLIT_P_PROJ + // Place the split point at the midpoint, unless the reference encroaches upon the resulting edges. + // In that case, the split point is taken such that its projection onto ap or bp is p. + Point s = midpoint(va->point(), vb->point()); + + typename Geom_traits::Side_of_bounded_sphere_3 side_of_bounded_sphere = geom_traits().side_of_bounded_sphere_3_object(); + + if (side_of_bounded_sphere(va->point(), mid, ref) == ON_BOUNDED_SIDE) { + Object result = Gt().intersect_3_object()(Line(va->point(), vb->point()), Plane(ref, va->point() - ref)); + CGAL_triangulation_assertion(assign(s, result)); + } + else if (side_of_bounded_sphere(vb->point(), mid, ref) == ON_BOUNDED_SIDE) { + Object result = Gt().intersect_3_object()(Line(va->point(), vb->point()), Plane(ref, vb->point() - ref)); + CGAL_triangulation_assertion(assign(s, result)); + } + + return s; +#endif + + // Place the split point at the projection of the reference onto ab, unless that is too close to a or b, or it is not on the edge. + // In that case, the split point is the midpoint between a and b. + Line l(va->point(), vb->point()); + Point s = l.projection(ref); + + // The following will enforce a minimum edge length. + if (has_larger_distance_to_point(s, ref, va->point()) || has_larger_distance_to_point(s, ref, vb->point())) + s = midpoint(va->point(), vb->point()); + + // Note that for inexact construction kernels, this assertion may fail for both projections and midpoints. + //CGAL_triangulation_assertion(l.has_on(s)); + + // This method doesn't actually produce a correct bounding ball, because the Steiner points may be inserted closer to the endpoint than the reference point.. + + return s; + +#ifdef SPLIT_SI + // Si's point segment recovery. + typedef Gt::Sphere_3 Sphere; + Point s; + Object result; + if (c == Vertex_handle()) { + // Rule 1. + FT sd, r = squared_distance(a->point(), b->point()); + if (4*(sd = squared_distance(a->point(), p->point())) < r) { + result = Gt().intersect_3_object()(Sphere(a->point(), sd), Segment(a->point(), b->point())); + } + else if (4*(sd = squared_distance(b->point(), p->point())) < r) { + result = Gt().intersect_3_object()(Sphere(b->point(), sd), Segment(a->point(), b->point())); + } + else { + Vertex_handle v = insert(midpoint(a->point(), b->point()), a->cell()); + return v; + } + if (!assign(s, result)) return Vertex_handle(); + } + else { + // Rule 2. + FT r = squared_distance(c->point(), p->point()); + result = Gt().intersect_3_object()(Sphere(c->point(), r), Segment(a->point(), b->point())); + if (!assign(s, result)) return Vertex_handle(); + + FT sd = squared_distance(s, p->point()); + if (sd > squared_distance(s, b->point())) { + //Rule 3. + if (4*sd < squared_distance(s, a->point())) { + result = Gt().intersect_3_object()(Sphere(s, sd), Segment(a->point(), b->point())); + if (!assign(s, result)) + return Vertex_handle(); + } + else { + s = midpoint(a->point(), s); + } + } + } + + return s; +#endif +} + +template < class Gt, class Tds, class Itag > +void Conforming_Delaunay_triangulation_3:: +conform_segment(Vertex_handle va, Vertex_handle vb, Points_iterator encr_begin, Points_iterator encr_end) { + CGAL_triangulation_precondition(va != vb); + // This insertion method is an adjustment of the segment insertion + // in the paper "Constrained Delaunay tetrahedral mesh generation and refinement" by Si (2010). + // The points of the tetrahedra intersected by ab (encroaching points) are collected; + // from these points, the point making the largest circumsphere with ab is the reference point. + // The reference point indicates where ab is split and the subsegments each get a subset + // of the encroaching points and then the process is repeated. + + // Set the initial part of the segment that is already in the triangulation to conforming. + Vertex_handle vi; + Cell_handle c; + int li, lj; + while (includes_edge(va, vb, vi, c, li, lj)) { + set_conforming(c, li, lj, true); + va = vi; + if (va == vb) + return; + } + + // If the edge is not in the triangulation, there must be encroaching points. + // Collect the points encroaching upon and close to ab. + if (encr_begin == encr_end) { + Point_list encroaching; + bool intersect = collect_encroaching(va, vb, std::back_inserter(encroaching), vi); + + if (!intersect && encroaching.empty()) { + // There are occurrences where none of the vertices of the cells intersected by ab encroach upon ab. + // In this case, we insert the midpoint of ab. + vi = insert(midpoint(va->point(), vb->point()), va->cell()); + vi->steiner() = true; + intersect = true; + } + + // Insert ai, and if needed ib. + conform_segment(va, vi, encroaching.begin(), encroaching.end()); + if (vi != vb) { + conform_segment(vi, vb, encr_begin, encr_end); + } + return; + } + + // Find the reference point. + // This is the encroaching point closest to the segment. + FT sr, radius2 = 0; + Point ref; + for (Points_iterator it = encr_begin; it != encr_end; ++it) { + CGAL_triangulation_assertion(!collinear(va->point(), *it, vb->point())); + sr = squared_radius(va->point(), *it, vb->point()); + if (radius2 < sr) { + radius2 = sr; + ref = *it; + } + } + CGAL_triangulation_assertion(radius2 > 0); + + // Split the edge based on its reference point. + Point s = construct_split_point(va, vb, ref); + Vertex_handle vs = insert(s, va->cell()); + vs->steiner() = true; + + // Sort the encroaching points to encroaching upon as, sb, and only ab. + Points_iterator after_as = sort_encroaching(encr_begin, encr_end, va->point(), vs->point()); + Points_iterator after_sb = sort_encroaching(after_as, encr_end, vs->point(), vb->point()); + + // Insert the segments on both sides of s as constraints. + conform_segment(va, vs, encr_begin, after_as); + conform_segment(vs, vb, after_as, after_sb); +} + +template < class Gt, class Tds, class Itag > +void Conforming_Delaunay_triangulation_3:: +mark_segment(Vertex_handle va, Vertex_handle vb) { + Vertex_handle vi; + Cell_handle c; + int li, lj; + while (va != vb) { + // All segments must be in the triangulation. + CGAL_triangulation_assertion_code(bool ok = ) + includes_edge(va, vb, vi, c, li, lj); + CGAL_triangulation_assertion(ok); + + mark_edge(c, li, lj); + + va = vi; + } +} + +template < class Gt, class Tds, class Itag > +void Conforming_Delaunay_triangulation_3:: +insert_marked(Vertex_handle va, Vertex_handle vb) { + // Insert the segment and then mark all parts of the segment. + insert_conforming(va, vb); + mark_segment(va, vb); +} + +template < class Gt, class Tds, class Itag > +void typename Conforming_Delaunay_triangulation_3:: +set_conforming(Cell_handle c, int li, int lj, bool C) { + // The cell is already in the triangulation, so the edge is already Delaunay. + Vertex_handle vi = c->vertex(li), vj = c->vertex(lj); + switch (dimension()) { + case 3: { + // Constrain the edge in all its incident cells. + Cell_circulator it = incident_cells(c, li, lj), start(it); + do { + it->set_conforming(it->index(vi), it->index(vj), C); + ++it; + } while (it != start); + break; + } + case 2: { + // Constrain the edge in both cells that share it. + Cell_handle n = c->neighbor(3-li-lj); + n->set_conforming(n->index(vi), n->index(vj), C); + } + case 1: { + // Constrain the edge in the cell itself. + c->set_conforming(li, lj, C); + break; + } + } +} + +template < class Gt, class Tds, class Itag > +void typename Conforming_Delaunay_triangulation_3:: +mark_edge(Cell_handle c, int li, int lj) { + // The cell is already in the triangulation, so the edge is already Delaunay. + Vertex_handle vi = c->vertex(li), vj = c->vertex(lj); + switch (dimension()) { + case 3: { + // Constrain the edge in all its incident cells. + Cell_circulator it = incident_cells(c, li, lj), start(it); + do { + it->mark(it->index(vi), it->index(vj)); + it++; + } while (it != start); + break; + } + case 2: { + // Constrain the edge in both cells that share it. + Cell_handle n = c->neighbor(3-li-lj); + n->mark(n->index(vi), n->index(vj)); + } + case 1: { + // Constrain the edge in the cell itself. + c->mark(li, lj); + break; + } + } +} + +template < class Gt, class Tds, class Itag > +template < class InputIterator > +int typename Conforming_Delaunay_triangulation_3:: +insert_conforming_loop(InputIterator begin, InputIterator end, Cell_handle hint) { + // Insert the points in the loop. + std::list segments; + insert_loop_points(begin, end, std::back_inserter(segments), hint); + + // Insert the segments in the loop. + for (std::list::iterator it = segments.begin(); it != segments.end(); ++it) + insert_conforming(*it); + + return segments.size(); +} + +template < class Gt, class Tds, class Itag > +template < class InputIterator > +int typename Conforming_Delaunay_triangulation_3:: +insert_marked_loop(InputIterator begin, InputIterator end, Cell_handle hint) { + // Insert the points in the loop. + std::list segments; + insert_loop_points(begin, end, std::back_inserter(segments), hint); + + // Insert the segments in the loop. + for (std::list::iterator it = segments.begin(); it != segments.end(); ++it) + insert_marked(*it); + + return segments.size(); +} + +template < class Gt, class Tds, class Itag > +template < class InputIterator > +int typename Conforming_Delaunay_triangulation_3:: +remove_conforming_loop(InputIterator begin, InputIterator end, Cell_handle hint) { + // Insert the points in the loop. + std::list segments; + insert_loop_points(begin, end, std::back_inserter(segments), hint); + + // Insert the segments in the loop. + for (std::list::iterator it = segments.begin(); it != segments.end(); ++it) + remove_conforming(*it); + + return segments.size(); +} + +template < class Gt, class Tds, class Itag > +template < class InputIterator, class OutputIterator > +void typename Conforming_Delaunay_triangulation_3:: +insert_loop_points(InputIterator begin, InputIterator end, OutputIterator out, Cell_handle hint) { + if (begin == end) + return; + + // Insert the points in the loop. + Vertex_handle prev = insert(*begin, hint), first(prev); + hint = prev->cell(); + for (InputIterator it = ++begin; it != end; ++it) { + Vertex_handle cur = insert(*it, hint); + if (prev != cur) + *out++ = Bi_vertex(prev, cur); + hint = cur->cell(); + prev = cur; + } + if (first != prev) + *out++ = Bi_vertex(prev, first); +} + +template < class Gt, class Tds, class Itag > +void Conforming_Delaunay_triangulation_3:: +restore_conforming(Cell_handle c, int li) { + // Restore the conforming edges of facet i of c, based on its neighboring cell. + Cell_handle n = c->neighbor(li); + Vertex_handle v1, v2; Edge e; + for (int j = 0; j < dimension()+1; ++j) { + if (j != li) { + e = opposite_edge(c, li, (li+j)&3); + v1 = c->vertex(e.second); + v2 = c->vertex(e.third); + if (n->is_marked(n->index(v1), n->index(v2))) + e.first->mark(e.second, e.third); + else + e.first->set_conforming(e.second, e.third, n->is_conforming(n->index(v1), n->index(v2))); + } + } +} + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Vertex_handle +Conforming_Delaunay_triangulation_3:: +insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj) { + switch (dimension()) { + case 3: { + Conflict_tester_3 tester(p, this); + Vertex_handle v = insert_in_conflict(p, lt, c, li, lj, tester, conforming_visitor); + return v; + } // dim 3 + case 2: { + Conflict_tester_2 tester(p, this); + Vertex_handle v = insert_in_conflict(p, lt, c, li, lj, tester, conforming_visitor); + return v; + } // dim 2 + default : + // dimension <= 1 + // Do not use the generic insert. + return Tr::insert(p, c); + } +} + +template < class Gt, class Tds, class Itag > +void Conforming_Delaunay_triangulation_3:: +insert_conforming(const Point& a, const Point& b, Cell_handle hint) { + // Insert the two points and then conform the segment between them. + Vertex_handle va = insert(a, hint); + Vertex_handle vb = insert(b, va->cell()); + if (va != vb) + insert_conforming(va, vb); +} + +template < class Gt, class Tds, class Itag > +template < class InputIterator > +int Conforming_Delaunay_triangulation_3:: +insert_conforming(InputIterator begin, InputIterator end, Cell_handle hint) { + // Insert the points. + std::list segments; + for (InputIterator it = begin; it != end; ++it) { + Vertex_handle va = insert(it->first, hint); + hint = va->cell(); + if (it->first != it->second) { + Vertex_handle vb = insert(it->second, hint); + hint = vb->cell(); + segments.push_back(Bi_vertex(va, vb)); + } + } + + // Insert the segments + for (std::list::iterator it = segment.begin(); it != segments.end(); ++it) + insert_conforming(*it); + + return segment.size(); +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +remove(Vertex_handle v) { + // The vertex cannot be removed if it is incident to a conforming edge. + if (are_there_incident_conforming(v)) + return false; + + // Collect the opposite facets as seen from the outside. + std::vector cells; cells.reserve(32); + incident_cells(v, std::back_inserter(cells)); + std::vector boundary; boundary.reserve(cells.size()); + + for (std::vector::const_iterator it = cells.begin(); it != cells.end(); ++it) { + Cell_handle n = (*it)->neighbor((*it)->index(v)); + boundary.push_back(Facet(n, n->index(*it))); + } + cells.clear(); + + // Remove the vertex. +#ifdef CGAL_DELAUNAY_3_OLD_REMOVE + if (dimension() == 3 && !test_dim_down(v)) { + remove_3D_ear(v); + } else { +#endif + cDT tmp; + Vertex_remover remover(tmp); + Tr::remove(v, remover); +#ifdef CGAL_DELAUNAY_3_OLD_REMOVE + } +#endif + + // Reinsert any conforming edges on the boundary of the retriangulated region. + for (std::vector::iterator it = boundary.begin(); it != boundary.end(); ++it) { + restore_conforming(mirror_facet(*it)); + } + + CGAL_triangulation_expensive_postcondition(is_valid()); + return true; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +remove_conforming(Cell_handle c, int li, int lj) { + // Mark the edge as non-conforming + set_conforming(c, li, lj, false); + return true; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +remove_conforming(Vertex_handle va, Vertex_handle vb) { + CGAL_triangulation_precondition(va != vb); + // All the parts of the segment ab that are constrained and + // not incident to a constrained facet are marked unconstrained. + // The return is true if the whole segment ab is covered by unconstrained edges; + // otherwise, returns false; + + // Check if first part of the segment is an edge in the triangulation. + Vertex_handle vi; + Cell_handle c; + int li, lj; + bool complete = true; + if (includes_edge(va, vb, vi, c, li, lj)) { + // Remove the first part. + if (!remove_conforming(c, li, lj)) + complete = false; + if (vi != vb) { + // Remove the next part. + if (!remove_conforming(vi, vb)) + complete = false; + } + return complete; + } + + return false; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +are_there_incident_conforming(Vertex_handle v) const { + if (dimension() > 0) { + std::vector edges; + edges.reserve(32); + finite_incident_edges(v, std::back_inserter(edges)); + for (std::vector::iterator it = edges.begin(); it != edges.end(); ++it) + if (is_conforming(*it)) + return true; + } + + return false; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +includes_edge(Vertex_handle va, Vertex_handle vb, Vertex_handle& vc, Cell_handle& c, int& li, int& lj) const { + CGAL_triangulation_precondition(!is_infinite(va) && !is_infinite(vb)); + + // Returns true if the line segment ab contains an edge e + // incident to a, false otherwise. + // If true, vc becomes the vertex of e distinct from a, + // c is a cell incident to e where e = (c, li, lj). + Vertex_handle v; + std::vector edges; edges.reserve(64); + finite_incident_edges(va, std::back_inserter(edges)); + for (std::vector::const_iterator it = edges.begin(); it != edges.end(); ++it) { + // Find the other vertex of the edge. + v = it->first->vertex(it->second); + if (v == va) v = it->first->vertex(it->third); + + if (is_infinite(v)) + continue; + + // Check if this edge is (part of) ab. + if (v == vb) { + vc = vb; + c = it->first; + li = it->second; + lj = it->third; + return true; + } + else if (sign((va->point() - v->point()) * (vb->point() - v->point())) == NEGATIVE && + collinear(va->point(), v->point(), vb->point())/* && + collinear_are_ordered_along_line(va->point(), v->point(), vb->point())*/) { + vc = v; + c = it->first; + li = it->second; + lj = it->third; + return true; + } + } + return false; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +is_valid(bool verbose, int level) const { + if (!tds().is_valid(verbose,level)) { + if (verbose) + std::cerr << "invalid data structure" << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + if (infinite_vertex() == Vertex_handle()) { + if (verbose) + std::cerr << "no infinite vertex" << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + int inf; + Cell_handle n; + Vertex_handle v1, v2; + switch (dimension()) { + case 3: { + for (All_cells_iterator it = all_cells_begin(); it != all_cells_end(); ++it) { + for (int i = 0; i < 4; ++i) { + if (!it->neighbor(i)->has_neighbor(it)) { + if (verbose) + std::cerr << "inconsistent neighbors " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + if (it->has_vertex(infinite_vertex(), inf)) { + for (int i = 0; i < 4; ++i) { + if (i != inf) { + if (is_conforming(it, i, inf)) { + if (verbose) + std::cerr << "conforming infinite edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + else { + is_valid_finite(it); + for (int i = 0; i < 4; ++i) { + n = it->neighbor(i); + if (!is_infinite(n->vertex(n->index(it)))) { + if (side_of_sphere(it, n->vertex(n->index(it))->point()) == ON_BOUNDED_SIDE) { + if (verbose) + std::cerr << "non-empty sphere " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + } + for (Finite_edges_iterator it = finite_edges_begin(); it != finite_edges_end(); ++it) { + v1 = it->first->vertex(it->second); + v2 = it->first->vertex(it->third); + bool conforming = is_conforming(*it); + + Cell_circulator cit = incident_cells(it->first, it->second, it->third, it->first), start(cit); + do { + if (cit->is_marked(cit->index(v1), cit->index(v2))) { + if (verbose) + std::cerr << "marked edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + if (cit->is_conforming(cit->index(v1), cit->index(v2)) != conforming) { + if (verbose) + std::cerr << "inconsistent conforming edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + cit++; + } while (cit != start); + } + break; + } + case 2: { + for (All_facets_iterator it = all_facets_begin(); it != all_facets_end(); ++it) { + for (int i = 0; i < 3; ++i) { + if (!it->first->neighbor(i)->has_neighbor(it->first)) { + if (verbose) + std::cerr << "inconsistent neighbors " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + if (it->first->has_vertex(infinite_vertex(), inf) && inf < 3) { + for (int i = 0; i < 3; ++i) { + if (i != inf) { + if (is_conforming(it->first, i, 3)) { + if (verbose) + std::cerr << "conforming infinite edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + else { + is_valid_finite(it->first); + for (int i = 0; i < 3; ++i) { + n = it->first->neighbor(i); + if(!is_infinite(n->vertex(n->index(it->first)))) { + if (side_of_circle(it->first, 3, n->vertex(n->index(it->first))->point()) == ON_BOUNDED_SIDE) { + if (verbose) + std::cerr << "non-empty circle " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + + Edge o = opposite_edge(it->first, i, 3); + if (o.first->is_marked(o.second, o.third)) { + if (verbose) + std::cerr << "marked edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + if (is_conforming(o) != + is_conforming(n, n->index(o.first->vertex(o.second)), n->index(o.first->vertex(o.third)))) { + if (verbose) + std::cerr << "inconsistent conforming edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + } + } + break; + } + case 1: { + for (Finite_edges_iterator it = finite_edges_begin(); it != finite_edges_end(); ++it) + is_valid_finite(it->first); + break; + } + } + if (verbose) + std::cerr << "valid conforming Delaunay triangulation" << std::endl; + return true; +} + +template < class Gt, class Tds, class Itag > +bool Conforming_Delaunay_triangulation_3:: +is_valid(Cell_handle c, bool verbose, int level) const { + if (!DT::is_valid(c, verbose, level)) { + CGAL_triangulation_assertion(false); + return false; + } + + switch (dimension()) { + case 3: { + for (int i = 0; i < 4; ++i) { + for (int j = i+1; j < 4; ++j) { + Vertex_handle v1 = c->vertex(i), v2 = c->vertex(j); + if (is_infinite(v1) || is_infinite(v2)) { + if (is_conforming(c, i, j)) { + if (verbose) + std::cerr << "conforming infinite edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + else { + bool conforming = is_conforming(c, i, j); + + Cell_circulator cit = incident_cells(c, i, j, c), start(cit); + do { + if (cit->is_marked(cit->index(v1), cit->index(v2))) { + if (verbose) + std::cerr << "marked edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + if (cit->is_conforming(cit->index(v1), cit->index(v2)) != conforming) { + if (verbose) + std::cerr << "inconsistent conforming edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + cit++; + } while (cit != start); + } + } + } + break; + } + case 2: { + for (int v=0; v<3; ++v) { + int i = (v+1)%3, + j = (v+2)%3; + + Vertex_handle v1 = c->vertex(i), v2 = c->vertex(j); + if (is_infinite(v1) || is_infinite(v2)) { + if (is_conforming(c, i, j)) { + if (verbose) + std::cerr << "conforming infinite edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + else { + if (c->is_marked(i, j)) { + if (verbose) + std::cerr << "marked edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + + Cell_handle n = c->neighbor(v); + if (c->is_conforming(i, j) != n->is_conforming(n->index(v1), n->index(v2))) { + if (verbose) + std::cerr << "inconsistent conforming edge " << std::endl; + CGAL_triangulation_assertion(false); + return false; + } + } + } + break; + } + } + if (verbose) + std::cerr << "valid Delaunay cell" << std::endl; + return true; +} + +template < class Gt, class Tds, class Itag > +template +bool Conforming_Delaunay_triangulation_3:: +collect_encroaching(Vertex_handle va, Vertex_handle vb, OutputIterator encroaching, Vertex_handle& vc) { + CGAL_triangulation_precondition(va != vb); + CGAL_triangulation_precondition(!is_infinite(va) && !is_infinite(vb)); + CGAL_triangulation_precondition(dimension() >= 2); + + // Fill encroaching with the points that are encroaching upon the segment ab. + // vc is set to the first vertex on ab starting from a; + // if ab intersects an edge or facet the intersection is inserted and true returned. + + // This is done by traversing cells from a towards b until either b is found, + // or any other collinear vertex, or an intersecting constraint. + Locate_type lt; int li, lj; + Encroaching_collecter cit(va, vb->point(), this, encroaching); + while (cit.has_next()) { + // Walk to the next cell. + ++cit; + + cit.traversed(lt, li, lj); + + // Stop when either a new vertex is reached, or a barrier (a conforming edge) is encountered. + if (lt == VERTEX) { + vc = cit->vertex(li); + if (va == vc) + continue; + return false; + } + + if (cit.barrier_hit()) { + vc = intersect(cit, lt, li, lj, va, vb); + return true; + } + } + + // The cell must contain the target. + CGAL_triangulation_assertion(cit->has_vertex(vb)); + vc = vb; + return false; +} + +template < class Gt, class Tds, class Itag > +typename Conforming_Delaunay_triangulation_3::Points_iterator +Conforming_Delaunay_triangulation_3:: +sort_encroaching(Points_iterator begin, Points_iterator end, const Point& a, const Point& b) const { + // Sorts the points encroaching upon ab to the front of the collection. + // Returns an iterator to the first point not encroaching upon ab. + Point temp; + Points_iterator middle = end; + for (Points_iterator it = begin; it != middle;) { + if (collinear(a, b, *it) || !is_encroaching(a, b, *it)) { + if (it != --middle) { + // Swap the current point to the back. + temp = *middle; + *middle = *it; + *it = temp; + } + } + else + ++it; + } + return middle; +} + +CGAL_END_NAMESPACE + +#endif // CGAL_CONFORMING_DELAUNAY_TRIANGULATION_3_H \ No newline at end of file diff --git a/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_cell_base_3.h b/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_cell_base_3.h new file mode 100644 index 00000000000..435baa885fe --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_cell_base_3.h @@ -0,0 +1,179 @@ +//A tetrahedral cell in a conforming Delaunay triangulation. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +// Based on CGAL/Triangulation_cell_base_3.h +// + +// Note that a constrained 3D Delaunay triangulation is partially conforming and guided. +// 1-dimensional constraints are specified by segments. +// 2-demensional constraints are specified by closure complete polygons. +// The constrained 3D DT must contain a collection of edges and faces that exactly cover +// each constrained segment or polygon. For example, a constrained edge will not necessarily +// be maintained as vertices are inserted, but after updating the triangulation to the vertex +// insertion, the constrained edge must be the union of one or more constrained edges. + +#ifndef CGAL_CONFORMING_TRIANGULATION_CELL_BASE_3_H +#define CGAL_CONFORMING_TRIANGULATION_CELL_BASE_3_H + +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + +template < typename Gt, typename Cb = Triangulation_cell_base_3 > +class Conforming_triangulation_cell_base_3: public Cb { +public: + typedef typename Cb::Vertex_handle Vertex_handle; + typedef typename Cb::Cell_handle Cell_handle; + + typedef Gt Geom_traits; + typedef typename Geom_traits::Point_3 Point; + + typedef Point* Point_container; + typedef Point* Point_iterator; + typedef const Point* Point_const_iterator; + + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Cb::template Rebind_TDS::Other Cb2; + typedef Conforming_triangulation_cell_base_3 Other; + }; // Rebind_TDS + +protected: + // The edge states. + bool EF[6]; + bool mm[6]; + +public: + Conforming_triangulation_cell_base_3(): Cb() {clear_conforming();} + Conforming_triangulation_cell_base_3(const Vertex_handle& v0, const Vertex_handle& v1, + const Vertex_handle& v2, const Vertex_handle& v3) + : Cb(v0, v1, v2, v3) {clear_conforming();} + Conforming_triangulation_cell_base_3(const Vertex_handle& v0, const Vertex_handle& v1, + const Vertex_handle& v2, const Vertex_handle& v3, + const Cell_handle& n0, const Cell_handle& n1, + const Cell_handle& n2, const Cell_handle& n3) + : Cb(v0, v1, v2, v3, n0, n1, n2, n3) {clear_conforming();} + + bool is_conforming(int i, int j) const {return EF[getEdgeIndex(i,j)];} + bool is_marked(int i, int j) const {return mm[getEdgeIndex(i,j)];} + + void set_conforming(bool c0, bool c1, bool c2, bool c3, bool c4, bool c5) { + EF[0] = c0; + EF[1] = c1; + EF[2] = c2; + EF[3] = c3; + EF[4] = c4; + EF[5] = c5; + clear_marked(); + } + void set_conforming(int i, int j, bool c) { + int index = getEdgeIndex(i,j); + EF[index] = c; + mm[index] = false; + } + void mark(int i, int j) { + int index = getEdgeIndex(i,j); + EF[index] = true; + mm[index] = true; + } + + void clear_conforming() {set_conforming(false, false, false, false, false, false);} + void clear_marked() { + mm[0] = false; + mm[1] = false; + mm[2] = false; + mm[3] = false; + mm[4] = false; + mm[5] = false; + } + virtual void clear() {clear_conforming();} + + bool has_conforming() const {return EF[0] || EF[1] || EF[2] || EF[3] || EF[4] || EF[5];} + + virtual void reorient() { + bool tmp = EF[1]; + EF[1] = EF[3]; + EF[3] = tmp; + tmp = mm[1]; + mm[1] = mm[3]; + mm[3] = tmp; + } + + virtual std::istream& read_cell(std::istream& is); + virtual std::ostream& write_cell(std::ostream& os) const; + +private: + int getEdgeIndex(int i, int j) const { + CGAL_triangulation_precondition(i>=0 && i<=3); + CGAL_triangulation_precondition(j>=0 && j<=3); + CGAL_triangulation_precondition(i != j); + return (i==0 || j==0) ? i+j-1 : i+j; + } +}; // Conforming_triangulation_cell_base_3 + +template < class Gt, class Cb > +inline std::istream& operator>>(std::istream& is, Conforming_triangulation_cell_base_3& c) { + is >> static_cast(c); + return c.read_cell(is); +} + +template < class Gt, class Cb > +inline std::ostream& operator<<(std::ostream &os, const Conforming_triangulation_cell_base_3& c) { + os << static_cast(c); + return c.write_cell(os); +} + +template < class Gt, class Cb > +std::istream& Conforming_triangulation_cell_base_3::read_cell(std::istream& is) { + char s; + for (int i = 0; i < 4; ++i) { + for (int j = i+1; j < 4; ++j) { + if (is_ascii(is)) + is >> s; + else + read(is, s); + if (s == 'C') + set_conforming(i, j, true); + } + } + + return is; +} + +template < class Gt, class Cb > +std::ostream& Conforming_triangulation_cell_base_3::write_cell(std::ostream& os) const { + for (int i = 0; i < 4; ++i) { + for (int j = i+1; j < 4; ++j) { + if (is_conforming(i, j)) { + os << "C"; + } + else + os << "N"; + if (is_ascii(os)) + os << ' '; + } + } + + return os; +} + +CGAL_END_NAMESPACE + +#endif // CGAL_CONFORMING_TRIANGULATION_CELL_BASE_3_H \ No newline at end of file diff --git a/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_vertex_base_3.h b/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_vertex_base_3.h new file mode 100644 index 00000000000..bd1902b9fed --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Conforming_triangulation_vertex_base_3.h @@ -0,0 +1,69 @@ +//A vertex in the conforming Delaunay triangulation structure. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +#ifndef CGAL_CONFORMING_TRIANGULATION_VERTEX_BASE_3_H +#define CGAL_CONFORMING_TRIANGULATION_VERTEX_BASE_3_H + +#include + +CGAL_BEGIN_NAMESPACE + +template < typename GT, typename Vb = Triangulation_vertex_base_3 > +class Conforming_triangulation_vertex_base_3: public Vb { + bool _s; + +public: + typedef typename Vb::Cell_handle Cell_handle; + typedef typename Vb::Point Point; + + template < typename TDS2 > struct Rebind_TDS { + typedef typename Vb::template Rebind_TDS::Other Vb2; + typedef Conforming_triangulation_vertex_base_3 Other; + }; + + Conforming_triangulation_vertex_base_3(): Vb(), _s(false) {} + Conforming_triangulation_vertex_base_3(const Point& p): Vb(p), _s(false) {} + Conforming_triangulation_vertex_base_3(const Point& p, bool s): Vb(p), _s(s) {} + Conforming_triangulation_vertex_base_3(const Point& p, Cell_handle c): Vb(p, c), _s(false) {} + Conforming_triangulation_vertex_base_3(const Point& p, bool s, Cell_handle c): Vb(p, c), _s(s) {} + Conforming_triangulation_vertex_base_3(Cell_handle c): Vb(c), _s(false) {} + + const bool& steiner() const {return _s;} + bool& steiner() {return _s;} +}; + +template < class GT, class Vb > +std::istream& +operator>>(std::istream &is, Conforming_triangulation_vertex_base_3 &v) { + is >> static_cast(v); + if (is_ascii(is)) is >> v.steiner(); + else read(is, v.steiner()); + return is; +} + +template < class GT, class Vb > +std::ostream& +operator<<(std::ostream &os, const Conforming_triangulation_vertex_base_3 &v) { + os << static_cast(v); + if (is_ascii(os)) os << ' '; + return os << v.steiner(); +} + +CGAL_END_NAMESPACE + +#endif // CGAL_CONFORMING_TRIANGULATION_VERTEX_BASE_3_H diff --git a/Conforming_triangulation_3/include/CGAL/Delaunay_triangulation_utils_3.h b/Conforming_triangulation_3/include/CGAL/Delaunay_triangulation_utils_3.h new file mode 100644 index 00000000000..c9442b475dc --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Delaunay_triangulation_utils_3.h @@ -0,0 +1,402 @@ +//Some additional utilities for the Delaunay triangulation structure. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +// Based on CGAL/Delaunay_triangulation_3.h +// +// For the love of something, why is the 3D triangulation package so badly written? + +#ifndef CGAL_DELAUNAY_TRIANGULATION_3_UTILS_H +#define CGAL_DELAUNAY_TRIANGULATION_3_UTILS_H + +#include + +#include +#include + +#include +#include "CCDT/Triangulation_segment_traverser_3.h" + +CGAL_BEGIN_NAMESPACE + +template < class DT > class Natural_neighbors_3; + +template < class Gt, class Tds > class Triangulation_cell_traverser_3; + +template < class Gt, + class Tds = Triangulation_data_structure_3 < Triangulation_vertex_base_3, Triangulation_cell_base_3 > > +class Delaunay_triangulation_utils_3: public Delaunay_triangulation_3 { + typedef Triangulation_data_structure Tds; + + typedef Delaunay_triangulation_utils_3 Self; + typedef Delaunay_triangulation_3 DT; + typedef Triangulation_3 Tr; + + friend class Natural_neighbors_3; + +public: + typedef Tds Triangulation_data_structure; + typedef Gt Geom_traits; + + typedef typename Gt::FT FT; + + typedef typename Gt::Point_3 Point; + typedef typename Gt::Vector_3 Vector; + typedef typename Gt::Segment_3 Segment; + typedef typename Gt::Line_3 Line; + typedef typename Gt::Triangle_3 Triangle; + typedef typename Gt::Tetrahedron_3 Tetrahedron; + typedef typename Gt::Plane_3 Plane; + + typedef typename DT::Cell_handle Cell_handle; + typedef typename DT::Vertex_handle Vertex_handle; + + typedef typename DT::Cell Cell; + typedef typename DT::Vertex Vertex; + typedef typename DT::Facet Facet; + typedef typename DT::Edge Edge; + + typedef std::pair Bi_vertex; + typedef Triple Tri_vertex; + + typedef std::pair Constraint_1; + typedef std::vector Constraint_2; + + typedef typename DT::Cell_circulator Cell_circulator; + typedef typename DT::Facet_circulator Facet_circulator; + typedef typename DT::Cell_iterator Cell_iterator; + typedef typename DT::Facet_iterator Facet_iterator; + typedef typename DT::Edge_iterator Edge_iterator; + typedef typename DT::Vertex_iterator Vertex_iterator; + + typedef typename DT::Finite_vertices_iterator Finite_vertices_iterator; + typedef typename DT::Finite_cells_iterator Finite_cells_iterator; + typedef typename DT::Finite_facets_iterator Finite_facets_iterator; + typedef typename DT::Finite_edges_iterator Finite_edges_iterator; + + typedef typename DT::All_cells_iterator All_cells_iterator; + + typedef Triangulation_segment_traverser_3 Cell_traverser; + + typedef typename DT::Locate_type Locate_type; + +#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_2 + using DT::cw; + using DT::ccw; + using DT::coplanar; + using DT::coplanar_orientation; + using DT::dimension; + using DT::finite_facets_begin; + using DT::finite_facets_end; + using DT::finite_vertices_begin; + using DT::finite_vertices_end; + using DT::finite_cells_begin; + using DT::finite_cells_end; + using DT::finite_edges_begin; + using DT::finite_edges_end; + using DT::geom_traits; + using DT::infinite_vertex; + using DT::is_valid; + using DT::mirror_vertex; + using DT::next_around_edge; + using DT::number_of_vertices; + using DT::orientation; + using DT::remove; + using DT::coplanar_side_of_bounded_circle; + using DT::side_of_oriented_sphere; + using DT::side_of_segment; + using DT::tds; + using DT::vertex_triple_index; +#endif + +protected: + typedef std::list Tri_vertex_collection; + typedef std::list Bi_vertex_collection; + +protected: + // Test whether a newly inserted point conflicts with the existing cells. + class Conflict_tester_3 { + protected: + const Point& _p; + const Self* _tr; + + public: + Conflict_tester_3(const Point& pt, const Self* tr): _p(pt), _tr(tr) {} + bool operator()(const Cell_handle c) const {return _tr->side_of_sphere(c, _p, true) == ON_BOUNDED_SIDE;} + Oriented_side compare_weight(const Point&, const Point&) const {return ZERO;} + bool test_initial_cell(Cell_handle) const {return true;} + }; // class Conflict_tester_3 + + // Test whether a newly inserted point conflicts with the existing cells. + class Conflict_tester_2 { + protected: + const Point& _p; + const Self* _tr; + + public: + Conflict_tester_2(const Point& pt, const Self* tr): _p(pt), _tr(tr) {} + bool operator()(const Cell_handle c) const {return _tr->side_of_circle(c, 3, _p, true) == ON_BOUNDED_SIDE;} + Oriented_side compare_weight(const Point&, const Point&) const {return ZERO;} + bool test_initial_cell(Cell_handle) const {return true;} + }; // class Conflict_tester_2 + + class Hidden_point_visitor { + public: + Hidden_point_visitor() {} + + template < class InputIterator > + void process_cells_in_conflict(InputIterator, InputIterator) const {} + void reinsert_vertices(Vertex_handle) {} + Vertex_handle replace_vertex(Cell_handle c, int index, const Point&) {return c->vertex(index);} + void hide_point(Cell_handle, const Point&) {} + }; // class Hidden_point_visitor + + Hidden_point_visitor hidden_point_visitor; + +public: + // Constructors for bi-vertices. + inline Bi_vertex to_bi_vertex(Cell_handle c, int li, int lj) const {return Bi_vertex(c->vertex(li), c->vertex(lj));} + inline Bi_vertex to_bi_vertex(const Edge& e) const {return to_bi_vertex(e.first, e.second, e.third);} + inline Bi_vertex sort_bi_vertex(Vertex_handle v1, Vertex_handle v2) const {if (v1 < v2) return Bi_vertex(v1, v2); return Bi_vertex(v2, v1);} + inline Bi_vertex sort_bi_vertex(const Bi_vertex& bv) const {return sort_bi_vertex(bv.first, bv.second);} + inline Bi_vertex sort_bi_vertex(Cell_handle c, int li, int lj) const {return sort_bi_vertex(c->vertex(li), c->vertex(lj));} + inline Bi_vertex sort_bi_vertex(const Edge& e) const {return sort_bi_vertex(e.first, e.second, e.third);} + + // Constructors for tri-vertices. + inline Tri_vertex to_tri_vertex(Cell_handle c, int li) const {return Tri_vertex(c->vertex((li+1)&3), c->vertex((li+2)&3), c->vertex((li+3)&3));} + inline Tri_vertex to_tri_vertex(const Facet& f) const {return to_tri_vertex(f.first, f.second);} + inline Tri_vertex sort_tri_vertex(Vertex_handle v1, Vertex_handle v2, Vertex_handle v3) const { + CGAL_triangulation_precondition(v1 != v2 && v1 != v3 && v2 != v3); + if (v1 < v2) { + if (v2 < v3) + return Tri_vertex(v3, v2, v1); + if (v3 < v1) + return Tri_vertex(v2, v1, v3); + return Tri_vertex(v2, v3, v1); + } + else { // v2 < v1 + if (v3 < v2) + return Tri_vertex(v1, v2, v3); + if (v1 < v3) + return Tri_vertex(v3, v1, v2); + return Tri_vertex(v1, v3, v2); + } + } + inline Tri_vertex sort_tri_vertex(const Tri_vertex& tv) const {return sort_tri_vertex(tv.first, tv.second, tv.third);} + inline Tri_vertex sort_tri_vertex(Cell_handle c, int li) const {return sort_tri_vertex(to_tri_vertex(c, li));} + inline Tri_vertex sort_tri_vertex(const Facet& f) const {return sort_tri_vertex(f.first, f.second);} + inline Tri_vertex oriented_tri_vertex(Cell_handle c, int li) const { + CGAL_triangulation_precondition(dimension() == 2 || dimension() == 3); + CGAL_triangulation_precondition((dimension() == 2 && li == 3) || (dimension() == 3 && li >= 0 && li <= 3) ); + CGAL_triangulation_precondition(!is_infinite(c, li)); + if ((li&1) == 0) + return Tri_vertex(c->vertex((li+2)&3), + c->vertex((li+1)&3), + c->vertex((li+3)&3)); + return Tri_vertex(c->vertex((li+1)&3), + c->vertex((li+2)&3), + c->vertex((li+3)&3)); + } + inline Tri_vertex oriented_tri_vertex(const Facet& f) const {return oriented_tri_vertex(f.first, f.second);} + +protected: + // tri_vertex methods similar to Triangulation_3. + Tri_vertex make_tri_vertex(const Facet& f) const { + return Tri_vertex(f.first->vertex(vertex_triple_index(f.second, 0)), + f.first->vertex(vertex_triple_index(f.second, 1)), + f.first->vertex(vertex_triple_index(f.second, 2))); + } + void make_canonical(Tri_vertex& t) const { + int i = (&*(t.first) < &*(t.second))? 0 : 1; + if (i == 0) + i = (&*(t.first) < &*(t.third))? 0 : 2; + else + i = (&*(t.second) < &*(t.third))? 1 : 2; + Vertex_handle tmp; + switch (i) { + case 0: return; + case 1: + tmp = t.first; + t.first = t.second; + t.second = t.third; + t.third = tmp; + return; + default: + tmp = t.first; + t.first = t.third; + t.third = t.second; + t.second = tmp; + } + } + + // PREDICATES + Orientation coplanar_orientation(const Point& p0, const Point& p1, const Point& p2, const Point& p) const { + return Gt().coplanar_orientation_3_object()(p0, p1, p2, p); + } + +public: + // CONSTRUCTORS + Delaunay_triangulation_utils_3(const Gt& gt = Gt()): DT(gt) {} + + // Copy constructor duplicates vertices and cells. + Delaunay_triangulation_utils_3(const Delaunay_triangulation_utils_3& tr): DT(tr) { + CGAL_triangulation_postcondition(is_valid()); + } + + // Create a 3D Delaunay triangulation from a number of points. + template < typename InputIterator > + Delaunay_triangulation_utils_3(InputIterator first, InputIterator last, const Gt& gt = Gt()): DT(gt) { + insert(first, last); + CGAL_triangulation_postcondition(is_valid()); + } + +public: + // INSERT POINT + + Vertex_handle insert(const Point& p, Cell_handle start = Cell_handle()); + virtual Vertex_handle insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj); + + template < class InputIterator > + int insert(InputIterator first, InputIterator last) { + int n = number_of_vertices(); + + std::vector points(first, last); + std::random_shuffle(points.begin(), points.end()); + spatial_sort(points.begin(), points.end(), geom_traits()); + + Cell_handle hint; + for (std::vector::const_iterator it = points.begin(), end = points.end(); it != end; ++it) + hint = insert(*it, hint)->cell(); + + return number_of_vertices() - n; + } + +public: + // MOVE + + Vertex_handle move_point(Vertex_handle v, const Point& p); + +public: + // Construct the edge of c opposite to ij. + Edge opposite_edge(Cell_handle c, int li, int lj) const { + CGAL_triangulation_precondition(li >= 0 && li < 4); + CGAL_triangulation_precondition(lj >= 0 && lj < 4); + CGAL_triangulation_precondition(li != lj); + + switch (6-li-lj) { // i + j + missing indices = 6. + case 1: return Edge(c, 0, 1); + case 2: return Edge(c, 0, 2); + case 3: return (li == 0 || lj == 0) + ? Edge(c, 1, 2) + : Edge(c, 0, 3); + case 4: return Edge(c, 1, 3); + case 5: return Edge(c, 2, 3); + } + + CGAL_triangulation_assertion(false); + return Edge(); + } + Edge opposite_edge(const Edge& e) const {return opposite_edge(e.first, e.second, e.third);} + + // Give the same facet as seen from the other side. + Facet mirror_facet(Cell_handle c, int li) const {return Facet(c->neighbor(li), c->neighbor(li)->index(c));} + Facet mirror_facet(const Facet& f) const {return mirror_facet(f.first, f.second);} + + // Construct a plane of a facet. + Plane plane(Cell_handle c, int li) const { + /* This should be removed to reduce orientation errors with the inexact kernel*/ + CGAL_triangulation_precondition(dimension() >= 2); + CGAL_triangulation_precondition(dimension() == 3 || li == 3); + CGAL_triangulation_precondition(li >= 0 && li <= 3); + CGAL_triangulation_precondition(!is_infinite(c, li)); + if ((li&1) == 0) + return Plane(c->vertex((li+2)&3)->point(), + c->vertex((li+1)&3)->point(), + c->vertex((li+3)&3)->point()); + return Plane(c->vertex((li+1)&3)->point(), + c->vertex((li+2)&3)->point(), + c->vertex((li+3)&3)->point()); + } + Plane plane(const Facet& f) const {return plane(f.first, f.second);} + + // Construct a line of an edge. + Line line(Cell_handle c, int li, int lj) const { + CGAL_triangulation_precondition(dimension() >= 1); + CGAL_triangulation_precondition(li >= 0 && li <= 3); + CGAL_triangulation_precondition(lj >= 0 && lj <= 3); + CGAL_triangulation_precondition(li != lj); + CGAL_triangulation_precondition(li + lj < dimension() * 2); + CGAL_triangulation_precondition(!is_infinite(c, li, lj)); + return Line(c->vertex(li)->point(), + c->vertex(lj)->point()); + } + Line line(const Edge& e) const {return line(e.first, e.second, e.third);} +}; // Delaunay_triangulation_utils_3 + +template < class Gt, class Tds > +typename Delaunay_triangulation_utils_3::Vertex_handle +Delaunay_triangulation_utils_3:: +insert(const Point& p, Cell_handle start) { + Locate_type lt; + int li, lj; + Cell_handle c = locate(p, lt, li, lj, start); + return insert(p, lt, c, li, lj); +} + +template < class Gt, class Tds > +typename Delaunay_triangulation_utils_3::Vertex_handle +Delaunay_triangulation_utils_3:: +insert(const Point& p, Locate_type lt, Cell_handle c, int li, int lj) { + switch (dimension()) { + case 3: { + Conflict_tester_3 tester(p, this); + Vertex_handle v = insert_in_conflict(p, lt, c, li, lj, tester, hidden_point_visitor); + return v; + }// dim 3 + case 2: { + Conflict_tester_2 tester(p, this); + return insert_in_conflict(p, lt, c, li, lj, tester, hidden_point_visitor); + }//dim 2 + default : + // dimension <= 1 + // Do not use the generic insert. + return Tr::insert(p, c); + } +} + +template < class Gt, class Tds > +typename Delaunay_triangulation_utils_3::Vertex_handle +Delaunay_triangulation_utils_3::move_point(Vertex_handle v, const Point& p) { + CGAL_triangulation_precondition(!is_infinite(v)); + CGAL_triangulation_expensive_precondition(is_vertex(v)); + + // Remember an incident vertex to restart the point location after the removal. + Cell_handle c = v->cell(); + Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0); + CGAL_triangulation_assertion(old_neighbor != v); + + if (!remove(v)) + return v; + + if (dimension() <= 0) + return insert(p); + return insert(p, old_neighbor->cell()); +} + +CGAL_END_NAMESPACE + +#endif // CGAL_DELAUNAY_TRIANGULATION_3_UTILS_H diff --git a/Conforming_triangulation_3/include/CGAL/Encroaching_collecter_3.h b/Conforming_triangulation_3/include/CGAL/Encroaching_collecter_3.h new file mode 100644 index 00000000000..a60becee1c8 --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Encroaching_collecter_3.h @@ -0,0 +1,123 @@ +//A traverser that collects points encroaching upon its path. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +// A class that traverses the cells of a triangulation by following a segment from source to target, +// while collecting points that encroach upon the segment. + +#ifndef CGAL_ENCROACHING_COLLECTER_3_H +#define CGAL_ENCROACHING_COLLECTER_3_H + +#include "Triangulation_segment_traverser_3.h" +//#include "Conforming_Delaunay_triangulation_3.h" + +CGAL_BEGIN_NAMESPACE + +template < class Gt, class Tds, class Itag > +class Conforming_Delaunay_triangulation_3; + +template < class Gt, class Tds, class Itag = No_intersection_tag, class Out = std::back_insert_iterator>>> +class Encroaching_collecter_3: public Triangulation_segment_traverser_3 { + typedef Encroaching_collecter_3 Self; + typedef Triangulation_segment_traverser_3 Base; + + typedef typename Gt::Plane_3 Plane; + +public: + typedef Conforming_Delaunay_triangulation_3 CDT; + + typedef typename Tds::Vertex Vertex; + typedef typename Tds::Cell Cell; + typedef typename Tds::Edge Edge; + typedef typename Tds::Facet Facet; + typedef typename Tds::Vertex_handle Vertex_handle; + typedef typename Tds::Cell_handle Cell_handle; + typedef typename Tds::Cell_iterator Cell_iterator; + typedef typename Tds::Cell_circulator Cell_circulator; + typedef typename Tds::Facet_iterator Facet_iterator; + typedef typename Tds::Facet_circulator Facet_circulator; + + typedef typename Gt::Point_3 Point; + typedef typename CDT::Locate_type Locate_type; + typedef typename CDT::Bi_vertex Bi_vertex; + +protected: + Out _out; + + Bi_vertex se; // These are used when the source lies on an edge. + +public: + Encroaching_collecter_3(): Base() {} + Encroaching_collecter_3(Vertex_handle v, const Point& t, const CDT* tr, Out output) + : Base(v, t, tr), _out(output) {if (_lt == CDT::EDGE) se = _tr->sort_bi_vertex(_pos, _li, _lj);} + Encroaching_collecter_3(const Point& s, const Point& t, const CDT* tr, Out output, Cell_handle hint = Cell_handle()) + : Base(s, t, tr, hint), _out(output) {if (_lt == CDT::EDGE) se = _tr->sort_bi_vertex(_pos, _li, _lj);} + + virtual bool barrier_hit() const {return _lt == CDT::EDGE && _pos->is_conforming(_li, _lj) && _tr->sort_bi_vertex(_pos, _li, _lj) != se;} + +protected: + // Traverse to the next cell along the line. + virtual void increment(); +}; + +template < class Gt, class Tds, class Itag, class Out > +void Encroaching_collecter_3::increment() { + // Walk to the next cell. + Base::increment(); + + // Check for encroaching points. + // Which points are checked is based on the type of simplex traversed. + const CDT* _cdt = dynamic_cast(_tr); + switch (_lt) { + case CDT::FACET: + for (int i = 0; i < 4; ++i) + if (i != _li && !_cdt->is_infinite(_pos->vertex(i)) && _cdt->is_encroaching(_source, _target, _pos->vertex(i)->point())) + *_out++ = _pos->vertex(i)->point(); + break; + case CDT::EDGE: { + Vertex_handle vi = _pos->vertex(_li), vj = _pos->vertex(_lj), vk; + Cell_handle c; + + // Check the vertices in the star around the edge. + Facet_circulator fit = _cdt->incident_facets(_pos, _li, _lj), start(fit); + do { + c = fit->first; + vk = c->vertex(6 - c->index(vi) - c->index(vj) - fit->second); + if (!_cdt->is_infinite(vk) && !collinear(_source, _target, vk->point()) && _cdt->is_encroaching(_source, _target, vk->point())) + *_out++ = vk->point(); + ++fit; + } while (fit != start); + + // Check the vertices of the edge. + if (_cdt->is_encroaching(_source, _target, vi->point())) + *_out++ = vi->point(); + if (_cdt->is_encroaching(_source, _target, vj->point())) + *_out++ = vj->point(); + break; + } + case CDT::VERTEX: + CGAL_triangulation_assertion(collinear(_source, _target, _pos->vertex(_li)->point())); + break; + default: + CGAL_triangulation_assertion(false); + break; + } +} + +CGAL_END_NAMESPACE + +#endif // CGAL_ENCROACHING_COLLECTER_3_H \ No newline at end of file diff --git a/Conforming_triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h b/Conforming_triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h new file mode 100644 index 00000000000..fb0506a89ef --- /dev/null +++ b/Conforming_triangulation_3/include/CGAL/Triangulation_segment_traverser_3.h @@ -0,0 +1,695 @@ +//A class that follows a straight line through a Delaunay triangulation structure. +//Copyright (C) 2012 Utrecht University +// +//This program is free software: 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. +// +//This program is distributed in the hope that it will be useful, +//but WITHOUT ANY WARRANTY; without even the implied warranty of +//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +//GNU General Public License for more details. +// +//You should have received a copy of the GNU General Public License +//along with this program. If not, see . +// +// Author(s): Thijs van Lankveld + +// A class that traverses the cells of a triangulation by following a segment from source to target. + +#ifndef CGAL_TRIANGULATION_SEGMENT_TRAVERSER_3_H +#define CGAL_TRIANGULATION_SEGMENT_TRAVERSER_3_H + +#include "Delaunay_triangulation_utils_3.h" + +#include + +CGAL_BEGIN_NAMESPACE + +template < class Gt, class Tds> +class Delaunay_triangulation_utils_3; + +template < class Gt, class Tds > +class Triangulation_segment_traverser_3: public Forward_circulator_base< typename Tds::Cell, std::ptrdiff_t, std::size_t > { + typedef Triangulation_segment_traverser_3 Self; + +public: + typedef Delaunay_triangulation_utils_3 Tr; + + typedef typename Tds::Vertex Vertex; + typedef typename Tds::Cell Cell; + typedef typename Tds::Edge Edge; + typedef typename Tds::Facet Facet; + typedef typename Tds::Vertex_handle Vertex_handle; + typedef typename Tds::Cell_handle Cell_handle; + typedef typename Tds::Cell_iterator Cell_iterator; + typedef typename Tds::Cell_circulator Cell_circulator; + typedef typename Tds::Facet_iterator Facet_iterator; + typedef typename Tds::Facet_circulator Facet_circulator; + + typedef typename Gt::Point_3 Point; + typedef typename Tr::Locate_type Locate_type; + +protected: + // The triangulation being traversed. + const Tr* _tr; + + // The source and target points of the traversal. + Point _source, _target; + + // The cell currently traversed and the previous one. + Cell_handle _pos, _prev; + + // Otherwise, they indicate the simplex through which this cell was entered, + // or the location of the source if it is in this cell. + int _li, _lj; + Locate_type _lt; + + // This bit signifies when a cell containing the target is found. + bool _done; + + // Where possible, facets are checked in a random order. + mutable Random rng; + +public: + Triangulation_segment_traverser_3(): _tr(NULL), _pos(), _prev(), _li(-1), _lj(-1), _lt(Tr::CELL), _done(true) {} + Triangulation_segment_traverser_3(Vertex_handle s, const Point& t, const Tr* tr); + Triangulation_segment_traverser_3(const Point& s, const Point& t, const Tr* tr, Cell_handle hint = Cell_handle()); + + Self& operator++(); + Self operator++(int); + Cell* operator->() {return &*_pos;} + Cell& operator*() {return *_pos;} + Cell_handle handle() {return _pos;} + operator const Cell_handle() const {return _pos;} + bool operator==(const Self& ct) const; + bool operator!=(const Self& ct) const; + + bool operator==(const Cell_handle& ch) const {return ch == _pos;} + bool operator!=(const Cell_handle& ch) const {return ch != _pos;} + + bool operator==(Nullptr_t CGAL_triangulation_assertion_code(n)) const; + bool operator!=(Nullptr_t n) const; + + // Check if a cell containing the target has been reached. + bool has_next() const {return !_done;} + + // Access the current and previous cells. + Cell_handle previous() const {return _prev;} + + // Traverse to a cell containing the target. + Cell_handle traverse(); + + // Get the type of simplex traversed last. + void traversed(Locate_type& lt, int& li, int& lj) const {lt = _lt; li = _li; lj = _lj;} + +protected: + // Traverse to the next cell along the line. + virtual void increment(); + +private: + inline int edgeIndex(int i, int j) const { + CGAL_triangulation_precondition(i>=0 && i<=3); + CGAL_triangulation_precondition(j>=0 && j<=3); + CGAL_triangulation_precondition(i != j); + return (i==0 || j==0) ? i+j-1 : i+j; + } + + inline Orientation coplanar_orientation(const Point& p, const Point& q, const Point& r, const Point& s) { + return Gt().coplanar_orientation_3_object()(p, q, r, s); + } +}; // class Triangulation_segment_traverser_3 + +template < class Gt, class Tds > +inline bool operator==(typename Tds::Cell_handle ch, Triangulation_segment_traverser_3 tci) {return tci == ch;} + +template < class Gt, class Tds > +inline bool operator!=(typename Tds::Cell_handle ch, Triangulation_segment_traverser_3 tci) {return tci != ch;} + +template < class Gt, class Tds > +Triangulation_segment_traverser_3::Triangulation_segment_traverser_3(Vertex_handle s, const Point& t, const Tr* tr) +: _tr(tr), _pos(), _prev(), _lj(-1), _lt(Tr::VERTEX), _done(false) { + CGAL_triangulation_precondition(!_tr->is_infinite(s)); + CGAL_triangulation_precondition(s->point() != t); + CGAL_triangulation_precondition(_tr->dimension() >= 2); + CGAL_triangulation_assertion(_tr->dimension() == 3 || _tr->plane(*_tr->finite_facets_begin()).has_on(_target)); + + _source = s->point(); + _target = t; + + _pos = s->cell(); + int inf; + if (_pos->has_vertex(_tr->infinite_vertex(), inf)) + _pos = _pos->neighbor(inf); + _li = _pos->index(s); + + CGAL_triangulation_postcondition(_pos != Cell_handle()); +} + +template < class Gt, class Tds > +Triangulation_segment_traverser_3::Triangulation_segment_traverser_3(const Point& s, const Point& t, const Tr* tr, Cell_handle hint) +: _tr(tr), _pos(), _prev(), _li(-1), _lj(-1), _done(false) { + CGAL_triangulation_precondition(s != t); + CGAL_triangulation_precondition(_tr->dimension() >= 2); + CGAL_triangulation_assertion(_tr->dimension() == 3 || _tr->plane(*_tr->finite_facets_begin()).has_on(_target)); + + _source = s; + _target = t; + + _pos = _tr->locate(s, _lt, _li, _lj, hint); + + CGAL_triangulation_postcondition(_pos != Cell_handle()); +} + +template < class Gt, class Tds > +inline Triangulation_segment_traverser_3& +Triangulation_segment_traverser_3::operator++() { + CGAL_triangulation_precondition(_pos != Cell_handle()); + increment(); + return *this; +} + +template < class Gt, class Tds > +inline Triangulation_segment_traverser_3 +Triangulation_segment_traverser_3::operator++(int) { + Self tmp(*this); + ++(*this); + return tmp; +} + +template < class Gt, class Tds > +inline bool Triangulation_segment_traverser_3::operator==(const Self& ct) const { + CGAL_triangulation_precondition(_pos != Cell_handle() && ct._pos != Cell_handle()); + return (_pos == ct._pos && + _prev == ct._prev && + _tr == ct._tr && + _li == ct._li && + _lj == ct._lj && + _lt == ct._lt && + _source == ct._source && + _target == ct._target && + _done == ct._done); +} + +template < class Gt, class Tds > +inline bool Triangulation_segment_traverser_3::operator!=(const Self& ct) const { + return !(*this == ct); +} + +template < class Gt, class Tds > +inline bool Triangulation_segment_traverser_3::operator==(Nullptr_t CGAL_triangulation_assertion_code(n)) const { + CGAL_triangulation_assertion(n == NULL); + return _pos == Cell_handle(); +} + +template < class Gt, class Tds > +inline bool Triangulation_segment_traverser_3::operator!=(Nullptr_t n) const { + CGAL_triangulation_assertion(n == NULL); + return !(*this == n); +} + +template < class Gt, class Tds > +inline typename Triangulation_segment_traverser_3::Cell_handle +Triangulation_segment_traverser_3::traverse() { + while (!_done) + increment(); + return _pos; +} + +template < class Gt, class Tds > +void Triangulation_segment_traverser_3::increment() { + CGAL_triangulation_precondition(!_done); + + // Walks to the next cell over a facet intersected by the line from source to target. + // This method is based on Triangulation_3::locate(). + switch (_tr->dimension()) { + case 3: { + // Infinite cells should be handled differently. + int inf; + if (_pos->has_vertex(_tr->infinite_vertex(), inf)) { + // If this cell was reached by traversal from a finite one, it must be the final cell. + Cell_handle fin = _pos->neighbor(inf); + if (fin == _prev) { + _done = true; + // Here, I do not change the _lt etc. because they were already set. + return; + } + + const Point* pts[4]; + for (int i = 0; i != 4; ++i) + if (i != inf) + pts[i] = &(_pos->vertex(i)->point()); + pts[inf] = &_target; + Orientation o[4]; + + // For the remembering stochastic walk, we start trying with a random index: + int li = rng.template get_bits<2>(); + + // Check if the target lies outside the convex hull. + if (orientation(*pts[0], *pts[1], *pts[2], *pts[3]) == POSITIVE) { + // The target lies in an infinite cell. + _done = true; + return; + + // Note that we do not traverse to other infinite cells. + } + + pts[inf] = &_source; + CGAL_triangulation_assertion(orientation(*pts[0], *pts[1], *pts[2], *pts[3]) == POSITIVE); + + // Check if the line enters an adjacent infinite cell. + // This occurs if the target lies on the other side of + // a plane through one of the finite edges and the source point. + for (int j = 0; j != 4; ++j, li = (li+1)&3) { + if (li == inf) { + o[li] = COPLANAR; + continue; + } + + Cell_handle next = _pos->neighbor(li); + if (next == _prev) { + o[li] = POSITIVE; + continue; + } + + const Point* backup = pts[li]; + pts[li] = &_target; + o[li] = orientation(*pts[0], *pts[1], *pts[2], *pts[3]); + + if (o[li] != NEGATIVE) { + pts[li] = backup; + continue; + } + + // The target lies behind the plane through the source and two finite vertices. + // Traverse to the incident infinite cell. + _prev = _pos; + _pos = next; + CGAL_triangulation_assertion(_tr->is_infinite(_pos)); + + _lt = Tr::FACET; + _li = _pos->index(_prev); + return; + } + + // The line enters the convex hull here (or lies on the finite facet). + _prev = _pos; + _pos = fin; + + // Check through which simplex the line traverses. + switch (o[0]+o[1]+o[2]+o[3]) { + case 3: + _lt = Tr::FACET; + _li = _pos->index(_prev); + return; + case 2: + _lt = Tr::EDGE; + for (int i = 0; i < 4; ++i) { + if (o[i] == COPLANAR && i != inf) { + Edge opp = _tr->opposite_edge(_prev, inf, i); + _li = _pos->index(_prev->vertex(opp.second)); + _lj = _pos->index(_prev->vertex(opp.third)); + return; + } + } + CGAL_triangulation_assertion(false); + return; + case 1: + _lt = Tr::VERTEX; + for (int i = 0; i < 4; ++i) { + if (o[i] == POSITIVE) { + _li = _pos->index(_prev->vertex(i)); + return; + } + } + CGAL_triangulation_assertion(false); + return; + default: + CGAL_triangulation_assertion(false); + return; + } + } + + const Point* pts[4] = {&(_pos->vertex(0)->point()), + &(_pos->vertex(1)->point()), + &(_pos->vertex(2)->point()), + &(_pos->vertex(3)->point())}; + + // We check in which direction the target lies + // by comparing its position relative to the planes through the + // source and the edges of the cell. + Orientation o[6]; + bool calc[6] = {false, false, false, false, false, false}; + + if (_lt == Tr::VERTEX) { + // The three planes through the vertex are set to coplanar. + for (int j = 0; j < 4; ++j) { + if (_li != j) { + int ij = edgeIndex(_li, j); + o[ij] = COPLANAR; + calc[ij] = true; + } + } + } + else if (_lt == Tr::EDGE) { + // The plane through the edge is set to coplanar. + int ij = edgeIndex(_li, _lj); + o[ij] = COPLANAR; + calc[ij] = true; + } + + // For the remembering stochastic walk, we start trying with a random facet: + int li = rng.template get_bits<2>(); + + CGAL_triangulation_assertion_code(bool incell = true;) + for (int k = 0; k < 4; ++k, li = (li+1)&3) { + Cell_handle next = _pos->neighbor(li); + if (next == _prev) + continue; + + // Check if the target is outside the cell. + const Point* backup = pts[li]; + pts[li] = &_target; + if (orientation(*pts[0], *pts[1], *pts[2], *pts[3]) != NEGATIVE) { + pts[li] = backup; + continue; + } + + CGAL_triangulation_assertion_code(incell = false;) + + // Check if the target is inside the pyramid. + int lj = rng.template get_bits<2>(); + int or = 0; + for (int l = 0; l < 4; ++l, lj = (lj+1)&3) { + if (li == lj) + continue; + + // We check the orientation of the target compared to the plane + // Through the source and the edge opposite of ij. + int oij = 5 - edgeIndex(li, lj); + if (!calc[oij]) { + const Point* backup2 = pts[lj]; + pts[lj] = &_source; + o[oij] = orientation(*pts[0], *pts[1], *pts[2], *pts[3]); + pts[lj] = backup2; + calc[oij] = true; + } + + or -= o[oij]; + if (o[oij] == POSITIVE) { + // The target is not inside the pyramid. + // Invert the planes. + // This can be safely done because either + // they were not calculated yet, + // or they will no longer be used. + for (int j = 0; j < 4; ++j) { + if (li == j) continue; + int oij = 5 - edgeIndex(li, j); + o[oij] = -o[oij]; + } + or = 0; + break; + } + } + + if (or == 0) { + // Either the target is not inside the pyramid, + // or the pyramid is degenerate. + pts[li] = backup; + continue; + } + + // The target is inside the pyramid. + _prev = _pos; + _pos = next; + + switch (or) { + case 3: + _lt = Tr::FACET; + _li = _pos->index(_prev); + return; + case 2: + _lt = Tr::EDGE; + for (int j = 0; j < 4; ++j) { + if (li != j && o[5-edgeIndex(li, j)] == COPLANAR) { + Edge opp = _tr->opposite_edge(_prev, li, j); + _li = _pos->index(_prev->vertex(opp.second)); + _lj = _pos->index(_prev->vertex(opp.third)); + return; + } + } + CGAL_triangulation_assertion(false); + return; + case 1: + _lt = Tr::VERTEX; + for (int j = 0; j < 4; ++j) { + if (li != j && o[5-edgeIndex(li, j)] == NEGATIVE) { + _li = _pos->index(_prev->vertex(j)); + return; + } + } + CGAL_triangulation_assertion(false); + return; + default: + CGAL_triangulation_assertion(false); + return; + } + } + + // The target lies inside this cell. + CGAL_triangulation_assertion(incell); + _done = true; + return; + } + case 2: { + int inf; + if (_pos->has_vertex(_tr->infinite_vertex(), inf) && inf < 3) { + // If this cell was reached by traversal from a finite one, it must be the final cell. + Cell_handle fin = _pos->neighbor(inf); + if (fin == _prev) { + _done = true; + return; + } + + // Check the neighboring cells. + if (coplanar_orientation(_source, _pos->vertex(_tr->ccw(inf))->point(), _pos->vertex(_tr->cw(inf))->point(), _target) == NEGATIVE) { + Cell_handle tmp = _pos->neighbor(_tr->cw(inf)); + _prev = _pos; + _pos = tmp; + return; + } + if (coplanar_orientation(_source, _pos->vertex(_tr->cw(inf))->point(), _pos->vertex(_tr->ccw(inf))->point(), _target) == NEGATIVE) { + Cell_handle tmp = _pos->neighbor(_tr->ccw(inf)); + _prev = _pos; + _pos = tmp; + return; + } + if (coplanar_orientation(_pos->vertex(_tr->ccw(inf))->point(), _pos->vertex(_tr->cw(inf))->point(), _source, _target) != POSITIVE) { + _prev = _pos; + _pos = fin; + return; + } + + // The target lies in this infinite cell. + _done = true; + return; + } + + const Point* pts[3] = {&(_pos->vertex(0)->point()), + &(_pos->vertex(1)->point()), + &(_pos->vertex(2)->point())}; + + switch (_lt) { + case Tr::VERTEX: { + // First we try the incident edges. + Orientation ocw = coplanar_orientation(*pts[_li], *pts[(_li+2)%3], *pts[(_li+1)%3], _target); + if (_pos->neighbor((_li+1)%3) != _prev && + ocw == NEGATIVE) { + Cell_handle tmp = _pos->neighbor((_li+1)%3); + _prev = _pos; + _pos = tmp; + _li = _pos->index(_prev->vertex(_li)); + return; + } + Orientation occw = coplanar_orientation(*pts[_li], *pts[(_li+1)%3], *pts[(_li+2)%3], _target); + if (_pos->neighbor((_li+2)%3) != _prev && + occw == NEGATIVE) { + Cell_handle tmp = _pos->neighbor((_li+2)%3); + _prev = _pos; + _pos = tmp; + _li = _pos->index(_prev->vertex(_li)); + return; + } + + // Then we try the opposite edge. + if (coplanar_orientation(*pts[(_li+1)%3], *pts[(_li+2)%3], *pts[_li], _target) == NEGATIVE) { + Cell_handle tmp = _pos->neighbor(_li); + _prev = _pos; + _pos = tmp; + + switch (ocw+occw) { + case 2: { + _lt = Tr::EDGE; + _li = _pos->index(_prev->vertex((_li+1)%3)); + _lj = _pos->index(_prev->vertex((_li+2)%3)); + return; + } + case 1: { + _lt = Tr::VERTEX; + if (ocw == COLLINEAR) _li = _pos->index(_prev->vertex((_li+2)%3)); + else _li = _pos->index(_prev->vertex((_li+1)%3)); + return; + } + default: + CGAL_triangulation_assertion(false); + return; + } + } + + // The target lies in this cell. + _done = true; + return; + } + case Tr::EDGE: { + int lk = 3-_li-_lj; + + if (_pos->neighbor(lk) != _prev) { + // Check the edge itself + switch (coplanar_orientation(*pts[_li], *pts[_lj], *pts[lk], _target)) { + case COLLINEAR: { + // The target lies in this cell. + _done = true; + return; + } + case NEGATIVE: { + // The target lies opposite of the edge. + Cell_handle tmp = _pos->neighbor(lk); + _prev = _pos; + _pos = tmp; + _li = _pos->index(_prev->vertex(_li)); + _lj = _pos->index(_prev->vertex(_lj)); + return; + } + default: + break; + } + } + + Orientation o = coplanar_orientation(_source, *pts[lk], *pts[_li], _target); + switch (o) { + case POSITIVE: { + // The ray passes through the edge ik. + if (coplanar_orientation(*pts[lk], *pts[_li], _source, _target) == NEGATIVE) { + Cell_handle tmp = _pos->neighbor(_lj); + _prev = _pos; + _pos = tmp; + + if (collinear(_source, *pts[_li], _target)) { + _lt = Tr::VERTEX; + _li = _pos->index(_prev->vertex(_li)); + } + else { + _lt = Tr::EDGE; + _li = _pos->index(_prev->vertex(_li)); + _lj = _pos->index(_prev->vertex(lk)); + } + return; + } + break; + } + default: { + // The ray passes through the edge jk. + if (coplanar_orientation(*pts[lk], *pts[_lj], _source, _target) == NEGATIVE) { + Cell_handle tmp = _pos->neighbor(_li); + _prev = _pos; + _pos = tmp; + + if (collinear(_source, *pts[_lj], _target)) { + _lt = Tr::VERTEX; + _li = _pos->index(_prev->vertex(_lj)); + } + else if (o == COLLINEAR) { + _lt = Tr::VERTEX; + _li = _pos->index(_prev->vertex(lk)); + } + else { + _lt = Tr::EDGE; + _li = _pos->index(_prev->vertex(lk)); + _lj = _pos->index(_prev->vertex(_lj)); + } + return; + } + break; + } + } + + // The target lies in this cell. + _done = true; + return; + } + case Tr::FACET: { + // We test its edges in a random order until we find a neighbor to go further + int li = rng.get_int(0, 3); + + Orientation o[3]; + bool calc[3] = {false, false, false}; + + for (int j = 0; j != 3; ++j, li = (li+1)%3) { + Cell_handle next = _pos->neighbor(li); + if (next == _prev) + continue; + + // The target should lie on the other side of the edge. + if (coplanar_orientation(*pts[(li+1)%3], *pts[(li+2)%3], *pts[li], _target) != NEGATIVE) + continue; + + // The target should lie inside the wedge. + if (!calc[(li+1)%3]) { + o[(li+1)%3] = coplanar_orientation(_source, *pts[(li+1)%3], *pts[(li+2)%3], _target); + calc[(li+1)%3] = true; + } + if (o[(li+1)%3] == NEGATIVE) continue; + + if (!calc[(li+2)%3]) { + o[(li+2)%3] = coplanar_orientation(_source, *pts[(li+2)%3], *pts[(li+1)%3], _target); + calc[(li+2)%3] = true; + } + if (o[(li+2)%3] == POSITIVE) continue; + + _prev = _pos; + _pos = next; + + switch (o[(li+1)%3]+o[(li+2)%3]) { + case 2: { + _lt = Tr::EDGE; + _li = _pos->index(_prev->vertex((li+1)%3)); + _lj = _pos->index(_prev->vertex((li+2)%3)); + return; + } + case 1: { + _lt = Tr::VERTEX; + if (o[(li+1)%3] == COLLINEAR) _li = _pos->index(_prev->vertex((li+1)%3)); + else _li = _pos->index(_prev->vertex((li+2)%3)); + return; + } + default: + CGAL_triangulation_assertion(false); + return; + } + } + + // The target lies in this cell. + _done = true; + return; + } + default: + CGAL_triangulation_assertion(false); + } + } + } +} + +CGAL_END_NAMESPACE + +#endif // CGAL_TRIANGULATION_SEGMENT_TRAVERSER_3_H \ No newline at end of file diff --git a/Conforming_triangulation_3/package_info/Conforming_triangulation_3/copyright b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/copyright new file mode 100644 index 00000000000..dea96f89cdd --- /dev/null +++ b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/copyright @@ -0,0 +1,2 @@ +Utrecht University (The Netherlands) & INRIA Sophia-Antipolis (France) + diff --git a/Conforming_triangulation_3/package_info/Conforming_triangulation_3/description.txt b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/description.txt new file mode 100644 index 00000000000..64a95593a98 --- /dev/null +++ b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/description.txt @@ -0,0 +1,3 @@ +Package Conforming_triangulation_3 : +conforming triangulation data structure +this can be wrapped around any CGAL::Triangulation_3, usually CGAL::Delaunay_triangulation_3 diff --git a/Conforming_triangulation_3/package_info/Conforming_triangulation_3/license.txt b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/license.txt new file mode 100644 index 00000000000..8bb8efcb72b --- /dev/null +++ b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) diff --git a/Conforming_triangulation_3/package_info/Conforming_triangulation_3/maintainer b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/maintainer new file mode 100644 index 00000000000..8dc617e8bb8 --- /dev/null +++ b/Conforming_triangulation_3/package_info/Conforming_triangulation_3/maintainer @@ -0,0 +1 @@ +Thijs van Lankveld