Merge branch 'Periodic_3_triangulation_3_inexact_locate-APelle'

- This small feature adds a method to locate point with inexact predicates in
  Periodic_3_triangulation_3.
  (cf. https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/Small_Features/Add_inexact_locate_in_Periodic_3_triangulation_3)
- Tested in CGAL-4.5-Ic-13
- Approved by the Release Manager.

Installation/changes.html

@@ -223,6 +223,12 @@ and <code>src/</code> directories).
     </li>
   </ul>
 
+  <h3>3D Periodic Triangulations</h3>
+  <ul>
+    <li>Add a method to locate point with inexact predicates.
+    </li>
+  </ul>
+
   <h3>3D Alpha Shapes</h3>
   <ul>
     <li> Add member functions in <code>CGAL::Alpha_shape_3</code> to give

Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/CGAL/Periodic_3_triangulation_3.h

@@ -785,6 +785,21 @@ search.
 Cell_handle 
 locate(const Point & query, Cell_handle start = Cell_handle()) const; 
 
+
+/*!
+Same as `locate()` but uses inexact predicates.
+
+This function returns a handle on a cell that is a good approximation of the exact
+location of `query`, while being faster. Note that it may return a handle on a cell
+whose interior does not contain `query`.
+
+Note that this function is available only if the %Cartesian coordinates of `query`
+are accessible with functions `x()`, `y()` and `z()`.
+*/
+Cell_handle
+inexact_locate(const Point & query, Cell_handle start = Cell_handle()) const;
+
+
 /*!
 The \f$ k\f$-face that contains `query` in its interior is 
 returned, by means of the cell returned together with `lt`, which 

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h

@@ -53,6 +53,11 @@
 #include <CGAL/internal/Exact_type_selector.h>
 #include <CGAL/NT_converter.h>
 
+#ifndef CGAL_NO_STRUCTURAL_FILTERING
+#include <CGAL/Triangulation_structural_filtering_traits.h>
+#include <CGAL/determinant.h>
+#endif // no CGAL_NO_STRUCTURAL_FILTERING
+
 namespace CGAL {
 
 template < class GT, class TDS > class Periodic_3_triangulation_3;
@@ -62,6 +67,28 @@ template < class GT, class TDS > std::istream& operator>>
 template < class GT, class TDS > std::ostream& operator<< 
     (std::ostream& os, const Periodic_3_triangulation_3<GT,TDS> &tr);
 
+#ifndef CGAL_NO_STRUCTURAL_FILTERING
+namespace internal {
+// structural filtering is performed only for EPIC
+struct Periodic_structural_filtering_3_tag {};
+struct No_periodic_structural_filtering_3_tag {};
+
+template <bool filter>
+struct Periodic_structural_filtering_selector_3 {
+#ifdef FORCE_STRUCTURAL_FILTERING
+  typedef Periodic_structural_filtering_3_tag  Tag;
+#else
+  typedef No_periodic_structural_filtering_3_tag  Tag;
+#endif
+};
+
+template <>
+struct Periodic_structural_filtering_selector_3<true> {
+  typedef Periodic_structural_filtering_3_tag  Tag;
+};
+}
+#endif // no CGAL_NO_STRUCTURAL_FILTERING
+
 /**\class Periodic_3_triangulation_3
  * 
  * \brief Implements functionality for computing in periodic space.
@@ -767,10 +794,120 @@ public:
   }
   //@}
 
+#ifdef CGAL_NO_STRUCTURAL_FILTERING
+  Cell_handle
+  periodic_locate(const Point & p, const Offset &o_p,
+	 Locate_type & lt, int & li, int & lj,
+	 Cell_handle start = Cell_handle()) const;
+#else // no CGAL_NO_STRUCTURAL_FILTERING
+#  ifndef CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS
+#    define CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS 2500
+#  endif // no CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS
+
+public:
+  Cell_handle
+  inexact_periodic_locate(const Point& p, const Offset &o_p,
+                 Cell_handle start = Cell_handle(),
+                 int max_num_cells = CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const;
+protected:
+  Cell_handle
+  exact_periodic_locate(const Point& p, const Offset &o_p,
+               Locate_type& lt,
+               int& li, int & lj,
+               Cell_handle start) const;
+
+  Cell_handle
+  generic_periodic_locate(const Point& p, const Offset &o_p,
+                 Locate_type& lt,
+                 int& li, int & lj,
+                 Cell_handle start,
+                 internal::Periodic_structural_filtering_3_tag) const {
+    return exact_periodic_locate(p, o_p, lt, li, lj, inexact_periodic_locate(p, o_p, start));
+  }
+
+  Cell_handle
+  generic_periodic_locate(const Point& p, const Offset &o_p,
+                 Locate_type& lt,
+                 int& li, int & lj,
+                 Cell_handle start,
+                 internal::No_periodic_structural_filtering_3_tag) const {
+    return exact_periodic_locate(p, o_p, lt, li, lj, start);
+  }
+
+  Orientation
+  inexact_orientation(const Point &p, const Point &q,
+                      const Point &r, const Point &s) const
+  {
+    const double px = to_double(p.x());
+    const double py = to_double(p.y());
+    const double pz = to_double(p.z());
+    const double qx = to_double(q.x());
+    const double qy = to_double(q.y());
+    const double qz = to_double(q.z());
+    const double rx = to_double(r.x());
+    const double ry = to_double(r.y());
+    const double rz = to_double(r.z());
+    const double sx = to_double(s.x());
+    const double sy = to_double(s.y());
+    const double sz = to_double(s.z());
+
+    const double pqx = qx - px;
+    const double pqy = qy - py;
+    const double pqz = qz - pz;
+    const double prx = rx - px;
+    const double pry = ry - py;
+    const double prz = rz - pz;
+    const double psx = sx - px;
+    const double psy = sy - py;
+    const double psz = sz - pz;
+
+    const double det = determinant(pqx, pqy, pqz,
+                                   prx, pry, prz,
+                                   psx, psy, psz);
+    if (det > 0) return POSITIVE;
+    if (det < 0) return NEGATIVE;
+    return ZERO;
+  }
+
+  Orientation
+  inexact_orientation(const Point &p, const Point &q,
+                      const Point &r, const Point &s,
+                      const Offset& o_p, const Offset& o_q,
+                      const Offset& o_r, const Offset& o_s) const
+  {
+    return inexact_orientation(construct_point(p, o_p),
+        construct_point(q, o_q),
+        construct_point(r, o_r),
+        construct_point(s, o_s));
+  }
+
+public:
+
+  Cell_handle
+  periodic_locate(const Point & p, const Offset &o_p,
+         Locate_type & lt, int & li, int & lj,
+         Cell_handle start = Cell_handle()) const
+  {
+    typedef Triangulation_structural_filtering_traits<Geometric_traits> TSFT;
+    typedef typename internal::Periodic_structural_filtering_selector_3<
+      TSFT::Use_structural_filtering_tag::value >::Tag Should_filter_tag;
+
+    return generic_periodic_locate(p, o_p, lt, li, lj, start, Should_filter_tag());
+  }
+
+  Cell_handle
+  inexact_locate(const Point& p,
+                 Cell_handle start = Cell_handle(),
+                 int max_num_cells = CGAL_PT3_STRUCTURAL_FILTERING_MAX_VISITED_CELLS) const
+  {
+	  return inexact_periodic_locate(p, Offset(), start, max_num_cells);
+  }
+#endif // no CGAL_NO_STRUCTURAL_FILTERING
+
 protected:
   /** @name Location helpers */ //@{
-  Cell_handle periodic_locate(const Point & p, const Offset &o_p,
+//  Cell_handle periodic_locate(const Point & p, const Offset &o_p,
-    Locate_type & lt, int & li, int & lj, Cell_handle start) const;
+//    Locate_type & lt, int & li, int & lj, Cell_handle start) const;
 
   Bounded_side side_of_cell(const Point & p, const Offset &off,
       Cell_handle c, Locate_type & lt, int & i, int & j) const;
@@ -1500,9 +1637,14 @@ periodic_triangle(const Cell_handle c, int i) const
   * returns a vertex (Cell_handle,li) if lt == VERTEX
   */
 template < class GT, class TDS >
-typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle
+inline typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle
-Periodic_3_triangulation_3<GT,TDS>::periodic_locate(
+Periodic_3_triangulation_3<GT,TDS>::
-    const Point & p, const Offset &o_p,
+#ifdef CGAL_NO_STRUCTURAL_FILTERING
+periodic_locate
+#else
+exact_periodic_locate
+#endif
+(const Point & p, const Offset &o_p,
     Locate_type & lt, int & li, int & lj, Cell_handle start) const {
   int cumm_off = 0;
   Offset off_query = o_p;
@@ -1672,6 +1814,130 @@ try_next_cell:
   return c;
 }
 
+
+#ifndef CGAL_NO_STRUCTURAL_FILTERING
+template < class GT, class TDS >
+typename Periodic_3_triangulation_3<GT,TDS>::Cell_handle
+Periodic_3_triangulation_3<GT,TDS>::
+inexact_periodic_locate(const Point& p, const Offset& o_p,
+               Cell_handle start,
+               int n_of_turns) const
+{
+	int cumm_off = 0;
+	Offset off_query = o_p;
+	if (number_of_vertices() == 0) {
+		return Cell_handle();
+	}
+	CGAL_triangulation_assertion(number_of_vertices() != 0);
+
+	if (start == Cell_handle()) {
+		start = cells_begin();
+	}
+
+	cumm_off = start->offset(0) | start->offset(1)
+	    		| start->offset(2) | start->offset(3);
+	if (is_1_cover() && cumm_off != 0) {
+		if (((cumm_off & 4) == 4) && (FT(2)*p.x()<(_domain.xmax()+_domain.xmin())))
+			off_query += Offset(1,0,0);
+		if (((cumm_off & 2) == 2) && (FT(2)*p.y()<(_domain.ymax()+_domain.ymin())))
+			off_query += Offset(0,1,0);
+		if (((cumm_off & 1) == 1) && (FT(2)*p.z()<(_domain.zmax()+_domain.zmin())))
+			off_query += Offset(0,0,1);
+	}
+
+	CGAL_triangulation_postcondition(start!=Cell_handle());
+	CGAL_triangulation_assertion(start->neighbor(0)->neighbor(
+			start->neighbor(0)->index(start))==start);
+	CGAL_triangulation_assertion(start->neighbor(1)->neighbor(
+			start->neighbor(1)->index(start))==start);
+	CGAL_triangulation_assertion(start->neighbor(2)->neighbor(
+			start->neighbor(2)->index(start))==start);
+	CGAL_triangulation_assertion(start->neighbor(3)->neighbor(
+			start->neighbor(3)->index(start))==start);
+
+	// We implement the remembering visibility/stochastic walk.
+
+	// Remembers the previous cell to avoid useless orientation tests.
+	Cell_handle previous = Cell_handle();
+	Cell_handle c = start;
+
+	// Now treat the cell c.
+try_next_cell:
+  --n_of_turns;
+	cumm_off =
+			c->offset(0) | c->offset(1) | c->offset(2) | c->offset(3);
+
+	bool simplicity_criterion = (cumm_off == 0) && (off_query.is_null());
+
+	// We know that the 4 vertices of c are positively oriented.
+	// So, in order to test if p is seen outside from one of c's facets,
+	// we just replace the corresponding point by p in the orientation
+	// test.  We do this using the arrays below.
+
+	Offset off[4];
+	const Point* pts[4] = { &(c->vertex(0)->point()),
+			&(c->vertex(1)->point()),
+			&(c->vertex(2)->point()),
+			&(c->vertex(3)->point()) };
+
+	if (!simplicity_criterion && is_1_cover() ) {
+		for (int i=0; i<4; i++) {
+			off[i] = int_to_off(c->offset(i));
+		}
+	}
+
+	if (!is_1_cover()) {
+		// Just fetch the vertices of c as points with offsets
+		for (int i=0; i<4; i++) {
+			pts[i] = &(c->vertex(i)->point());
+			off[i] = get_offset(c,i);
+		}
+	}
+
+	for (int i=0; i != 4; ++i) {
+		Cell_handle next = c->neighbor(i);
+		if (previous == next) {
+			continue;
+		}
+
+		// We temporarily put p at i's place in pts.
+		const Point* backup = pts[i];
+		pts[i] = &p;
+
+		if (simplicity_criterion && is_1_cover() ) {
+			if ( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3]) != NEGATIVE ) {
+				pts[i] = backup;
+				continue;
+			}
+		}
+		else {
+			Offset backup_off;
+
+			backup_off = off[i];
+			off[i] = off_query;
+
+			if ( inexact_orientation(*pts[0], *pts[1], *pts[2], *pts[3],
+					off[0], off[1], off[2], off[3]) != NEGATIVE ) {
+				pts[i] = backup;
+				off[i] = backup_off;
+				continue;
+			}
+		}
+
+		// Test whether we need to adapt the offset of the query point.
+		// This means, if we get out of the current cover.
+		off_query = combine_offsets(off_query, get_neighbor_offset(c,i,next));
+		previous = c;
+		c = next;
+		if (n_of_turns)
+		  goto try_next_cell;
+	}
+
+	return c;
+}
+#endif
+
+
 /**
  * returns
  * ON_BOUNDED_SIDE if p inside the cell

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/include/CGAL/_test_cls_periodic_3_triangulation_3.h

@@ -419,6 +419,11 @@ _test_cls_periodic_3_triangulation_3(const PeriodicTriangulation &,
   assert(!PT3.are_equal(Facet(ch,1),nb,i));
 
   std::cout<<"Point location"<<std::endl;
+  ch = PT3_deg.inexact_locate(Point(0.5,0.5,0.5));
+  assert( PT3_deg.tetrahedron(PT3_deg.periodic_tetrahedron(ch))
+      == PT3_deg.geom_traits().construct_tetrahedron_3_object()(
+    Point(0,0,0), Point(0,2,0), Point(0,0,2), Point(2,0,0)) );
+
   ch = PT3_deg.locate(Point(0.5,0.5,0.5));
   assert( PT3_deg.tetrahedron(PT3_deg.periodic_tetrahedron(ch))
       == PT3_deg.geom_traits().construct_tetrahedron_3_object()(