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