Fix predicate/construction inconsistencies in P3T3 canonicalization functions

There are two use cases:
- either we care about the shifted position (for example for P3M3 where
we will insert the shifted position)
- either we don't care about the shifted position, and we are ok with
point + shifting offset (for example in P3M3 predicates to determine
is_bad)

In the first case, we cannot determine the shift using predicates,
otherwise we could have an inconsitency in the final result:
the predicates say the point is in, but once constructed it is not.

So this commit distinguishes between both. When we care about the
actual shifted position, we construct the point. There might be
numerical errors if we are not using exact constructions, but it
does not really matter.

What should be done better:
- use compare_x/y/z_3 instead of <
- handle the case where the numerical errors are such that you
get a really silly point far from the truth. Maybe this should
all be done in EPECK. There is something like this in the
history of canonicalize_helper.h ....

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp

@@ -55,7 +55,7 @@ public:
 
   return_type operator()(const Point_3& p) const
   {
-    const Point_3 cp = P3T3::internal::robust_canonicalize_point(p, m_gt);
+    const Point_3 cp = P3T3::internal::construct_canonical_point(p, m_gt);
     CGAL::Bounded_side res = m_sotm(cp);
     return ((res == ON_BOUNDED_SIDE) ? 1 : 2); // set a region to '0' if it is not to be meshed
   }

Periodic_3_mesh_3/include/CGAL/Periodic_3_function_wrapper.h

@@ -41,7 +41,7 @@ public:
 
   /// Operator ()
   return_type operator()(const Point_3& p) const {
-    return f_(P3T3::internal::robust_canonicalize_point(p, gt_));
+    return f_(P3T3::internal::construct_canonical_point(p, gt_));
   }
 
 private:

Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h

@@ -193,14 +193,14 @@ public:
 
   Bare_point canonicalize_point(const Bare_point& p) const
   {
-    return P3T3::internal::robust_canonicalize_point(p, geom_traits());
+    return P3T3::internal::construct_canonical_point(p, geom_traits());
   }
 
   // @fixme it might be dangerous to call robust_canonicalize() without also changing
   // <p, offset> = construct_periodic_point(p) (lack of consistency in the result)
   Weighted_point canonicalize_point(const Weighted_point& p) const
   {
-    return P3T3::internal::robust_canonicalize_point(p, geom_traits());
+    return P3T3::internal::construct_canonical_point(p, geom_traits());
   }
 
   // 1-cover, so we can take a const&
@@ -430,9 +430,8 @@ public:
     typename Geom_traits::Compute_squared_distance_3 compute_sd =
       geom_traits().compute_squared_distance_3_object();
 
-    bool used_exact = false;
-    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits());
-    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, geom_traits());
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, geom_traits());
 
     Offset min_off;
 
@@ -474,8 +473,7 @@ public:
     typename Geom_traits::Construct_point_3 cp =
       geom_traits().construct_point_3_object();
 
-    bool used_exact = false;
-    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, geom_traits());
 
     Offset min_off;
     Offset null_offset(0,0,0);
@@ -550,10 +548,9 @@ public:
       geom_traits().compare_power_distance_3_object();
 
     // Compute the offsets that would bring p, q, and r into the canonical domain
-    bool used_exact = false;
-    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits());
-    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(cp(q), used_exact, geom_traits());
-    std::pair<Bare_point, Offset> pp_r = P3T3::internal::construct_periodic_point(cp(r), used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, geom_traits());
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(cp(q), geom_traits());
+    std::pair<Bare_point, Offset> pp_r = P3T3::internal::construct_periodic_point(cp(r), geom_traits());
 
     // To compare pp(p, q) to pp(p, r), we first need to know the best offsets that minimize these distances
     auto get_offset_minimizing_power_product = [&](const Weighted_point& wp,

Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h

@@ -765,22 +765,15 @@ public:
   }
 
   // same as the base construct_periodic_point(), but for weighted points
-  Periodic_weighted_point construct_periodic_weighted_point(const Weighted_point& p,
-                                                            bool& had_to_use_exact) const
+  Periodic_weighted_point construct_periodic_weighted_point(const Weighted_point& p) const
   {
     const Bare_point& bp = geom_traits().construct_point_3_object()(p);
-    const Periodic_bare_point pbp = Tr_Base::construct_periodic_point(bp, had_to_use_exact);
+    const Periodic_bare_point pbp = Tr_Base::construct_periodic_point(bp);
     return std::make_pair(geom_traits().construct_weighted_point_3_object()(
                             pbp.first, p.weight()),
                           pbp.second);
   }
 
-  Periodic_weighted_point construct_periodic_weighted_point(const Weighted_point& p) const
-  {
-    bool useless = false;
-    return construct_periodic_weighted_point(p, useless);
-  }
-
 public:
   /** @name Geometric access functions */
   // The following functions change the position of a vertex.

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h

@@ -870,18 +870,11 @@ public:
     return construct_point(pp.first, pp.second);
   }
 
-  Periodic_point_3 construct_periodic_point(const Point_3& p,
-                                            bool& had_to_use_exact) const
+  Periodic_point_3 construct_periodic_point(const Point_3& p) const
   {
     // The function is a different file to be able to be used where there is
     // no triangulation (namely, the domains of Periodic_3_mesh_3).
-    return ::CGAL::P3T3::internal::construct_periodic_point(p, had_to_use_exact, geom_traits());
-  }
-
-  Periodic_point_3 construct_periodic_point(const Point_3& p) const
-  {
-    bool useless = false;
-    return construct_periodic_point(p, useless);
+    return ::CGAL::P3T3::internal::construct_periodic_point(p, geom_traits());
   }
 
   // ---------------------------------------------------------------------------

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h

@@ -18,7 +18,7 @@
 // Although 'canoncalize_point()' is used by Periodic_3_mesh_3, this file is here
 // (in Periodic_3_triangulation_3) because of 'construct_periodic_point()',
 // which is a function used in P3T3.h and also needed by 'canonicalize_point()'.
-// However, P3M3 needs 'canoncalize_point()' without having access to a triangulation
+// However, P3M3 needs 'canonicalize_point()' without having access to a triangulation
 // and to avoid duplicating it, the function is here.
 
 // Geom_traits must be a model of the concept 'P3T3Traits' for 'construct_periodic_point()'.
@@ -44,7 +44,6 @@ namespace internal {
 template <typename Gt_>
 std::pair<typename Gt_::Point_3, typename Gt_::Periodic_3_offset_3>
 construct_periodic_point(const typename Gt_::Point_3& p,
-                         bool& encountered_issue,
                          const Gt_& gt)
 {
   typedef Gt_                                                  Geom_traits;
@@ -52,13 +51,13 @@ construct_periodic_point(const typename Gt_::Point_3& p,
   typedef typename Geom_traits::Periodic_3_offset_3            Offset;
   typedef typename Geom_traits::Iso_cuboid_3                   Iso_cuboid;
 
-  const Iso_cuboid& domain = gt.get_domain();
-
-  // Use these rather than Construct_point_3 to avoid construction inaccuracies
+  // this assumes exact predicates, otherwise we could be looping
   typename Geom_traits::Compare_x_3 cmp_x3 = gt.compare_x_3_object();
   typename Geom_traits::Compare_y_3 cmp_y3 = gt.compare_y_3_object();
   typename Geom_traits::Compare_z_3 cmp_z3 = gt.compare_z_3_object();
 
+  const Iso_cuboid& domain = gt.get_domain();
+
   // Check if p lies within the domain. If not, translate.
   if(!(p.x() < domain.xmin()) && p.x() < domain.xmax() &&
      !(p.y() < domain.ymin()) && p.y() < domain.ymax() &&
@@ -67,108 +66,45 @@ construct_periodic_point(const typename Gt_::Point_3& p,
     return std::make_pair(p, Offset());
   }
 
-  // Numerical approximations might create inconsistencies between the constructions
-  // and the comparisons. For example in a cubic domain of size 2:
-  // 1. initial point: P(2+1e-17, 0, 0)
-  // 2. the function computes an offset(1, 0, 0),
-  // 3. P + (-1, 0, 0) * domain_size constructs Q(-1e-17, 0, 0) // numerical approximation
-  // 4. the function computes an offset of (-1, 0, 0)
-  // 5. Q + (1, 0, 0) * domain_size constructs (2+1e-17, 0, 0) (that is P)
-  // And the function is looping...
-  //
-  // If this is happening the 'Last_change' enum will break this infinite
-  // loop and return the wrong point and the 'encountered_issue' bool will be
-  // set to 'true'. An exact version of this function is then be called.
-
-  enum Last_change
-  {
-    NO_LAST_CHANGE,
-    INCREASED_X, DECREASED_X, INCREASED_Y, DECREASED_Y, INCREASED_Z, DECREASED_Z
-  };
-
-  Last_change lc = NO_LAST_CHANGE;
-  bool in = false;
+  Point domain_m(domain.xmin(), domain.ymin(), domain.zmin());
+  Point domain_M(domain.xmax(), domain.ymax(), domain.zmax());
 
   Offset transl(0, 0, 0);
   const Offset null_off(0, 0, 0);
 
-  Point domain_m(domain.xmin(), domain.ymin(), domain.zmin());
-  Point domain_M(domain.xmax(), domain.ymax(), domain.zmax());
-
-  while(!in)
+  for(;;)
   {
     if(cmp_x3(p, domain_m, transl, null_off) == SMALLER)
-    {
-      if(lc == DECREASED_X) // stuck in a loop
-        break;
-
-      lc = INCREASED_X;
       transl.x() += 1;
-    }
     else if(cmp_y3(p, domain_m, transl, null_off) == SMALLER)
-    {
-      if(lc == DECREASED_Y) // stuck in a loop
-        break;
-
-      lc = INCREASED_Y;
       transl.y() += 1;
-    }
     else if(cmp_z3(p, domain_m, transl, null_off) == SMALLER)
-    {
-      if(lc == DECREASED_Z) // stuck in a loop
-        break;
-
-      lc = INCREASED_Z;
       transl.z() += 1;
-    }
     else if(!(cmp_x3(p, domain_M, transl, null_off) == SMALLER))
-    {
-      if(lc == INCREASED_X) // stuck in a loop
-        break;
-
-      lc = DECREASED_X;
       transl.x() -= 1;
-    }
     else if(!(cmp_y3(p, domain_M, transl, null_off) == SMALLER))
-    {
-      if(lc == INCREASED_Y) // stuck in a loop
-        break;
-
-      lc = DECREASED_Y;
       transl.y() -= 1;
-    }
     else if(!(cmp_z3(p, domain_M, transl, null_off) == SMALLER))
-    {
-      if(lc == INCREASED_Z) // stuck in a loop
-        break;
-
-      lc = DECREASED_Z;
       transl.z() -= 1;
-    }
     else
-    {
-      in = true;
-    }
+      break;
   }
 
   std::pair<Point, Offset> pp(p, transl);
 
-  if(cmp_x3(p, domain_m, transl, null_off) == SMALLER || // < min
-     cmp_y3(p, domain_m, transl, null_off) == SMALLER ||
-     cmp_z3(p, domain_m, transl, null_off) == SMALLER ||
-     !(cmp_x3(p, domain_M, transl, null_off) == SMALLER) || // >= max
-     !(cmp_y3(p, domain_M, transl, null_off) == SMALLER) ||
-     !(cmp_z3(p, domain_M, transl, null_off) == SMALLER))
-  {
-    encountered_issue = true;
-  }
+  CGAL_postcondition(!(cmp_x3(p, domain_m, transl, null_off) == SMALLER) &&
+                     !(cmp_y3(p, domain_m, transl, null_off) == SMALLER) &&
+                     !(cmp_z3(p, domain_m, transl, null_off) == SMALLER) &&
+                       cmp_x3(p, domain_M, transl, null_off) == SMALLER &&
+                       cmp_y3(p, domain_M, transl, null_off) == SMALLER &&
+                       cmp_z3(p, domain_M, transl, null_off) == SMALLER);
 
   return pp;
 }
 
 template <typename Gt_>
 typename Gt_::Point_3
-constrain_to_canonical_domain(const typename Gt_::Point_3& p,
+clamp_to_canonical_domain(const typename Gt_::Point_3& p,
                           const Gt_& gt)
 {
   typedef Gt_                                         Geom_traits;
@@ -196,31 +132,114 @@ constrain_to_canonical_domain(const typename Gt_::Point_3& p,
 /// instance of the same bare point that lives inside the base domain
 template <typename Gt_>
 typename Gt_::Point_3
-robust_canonicalize_point(const typename Gt_::Point_3& p,
+construct_canonical_point(const typename Gt_::Point_3& p,
                           const Gt_& gt)
 {
   typedef Gt_                                                  Geom_traits;
-  typedef typename Geom_traits::Point_3               Bare_point;
   typedef typename Geom_traits::Periodic_3_offset_3            Offset;
   typedef typename Geom_traits::Iso_cuboid_3                   Iso_cuboid;
 
+  const Iso_cuboid& domain = gt.get_domain();
+
+  // Check if p lies within the domain. If not, translate.
+  if(!(p.x() < domain.xmin()) && p.x() < domain.xmax() &&
+     !(p.y() < domain.ymin()) && p.y() < domain.ymax() &&
+     !(p.z() < domain.zmin()) && p.z() < domain.zmax())
+  {
+    return p;
+  }
+
   typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
 
-  const Iso_cuboid& domain = gt.get_domain();
-  if(p.x() >= domain.xmin() && p.x() < domain.xmax() &&
-     p.y() >= domain.ymin() && p.y() < domain.ymax() &&
-     p.z() >= domain.zmin() && p.z() < domain.zmax())
-    return p;
+  // Numerical approximations might create inconsistencies between the constructions
+  // and the comparisons. For example in a cubic domain of size 2:
+  // 1. initial point: P(2+1e-17, 0, 0)
+  // 2. the function computes an offset(1, 0, 0),
+  // 3. P + (-1, 0, 0) * domain_size constructs Q(-1e-17, 0, 0) // numerical approximation
+  // 4. the function computes an offset of (-1, 0, 0)
+  // 5. Q + (1, 0, 0) * domain_size constructs (2+1e-17, 0, 0) (that is P)
+  // And the function is looping...
+  //
+  // If this is happening the 'Last_change' enum will break this infinite loop,
+  // and the function will snap the point to the domain.
 
-  bool encountered_issue = false;
-  std::pair<Bare_point, Offset> pbp = construct_periodic_point(p, encountered_issue, gt);
-  Bare_point canonical_p = cp(pbp.first /*point*/, pbp.second /*offset*/);
-
-  if(encountered_issue)
+  enum Last_change
   {
-    // If we encountered an issue, there's no guarantee that the double construction gives a point
-    // in the domain (even if we computed it exactly beforehand). So, forcefully put it into the domain.
-    canonical_p = constrain_to_canonical_domain(canonical_p, gt);
+    NO_LAST_CHANGE,
+    INCREASED_X, DECREASED_X, INCREASED_Y, DECREASED_Y, INCREASED_Z, DECREASED_Z
+  };
+
+  Last_change lc = NO_LAST_CHANGE;
+
+  Offset transl(0, 0, 0);
+  Point_3 canonical_p = p;
+
+  for(;;)
+  {
+    // Here we actually construct because we are using a kernel with exact constructions
+    // so it is faster than using <Object, Offset> APIs which have to perform constructions
+    // within the predicates every single time
+    canonical_p = cp(p, transl);
+
+    if(canonical_p.x() < domain.xmin())
+    {
+      if(lc == DECREASED_X) // stuck in a loop
+        break;
+
+      lc = INCREASED_X;
+      transl.x() += 1;
+    }
+    else if(canonical_p.y() < domain.xmin())
+    {
+      if(lc == DECREASED_Y) // stuck in a loop
+        break;
+
+      lc = INCREASED_Y;
+      transl.y() += 1;
+    }
+    else if(canonical_p.z() < domain.xmin())
+    {
+      if(lc == DECREASED_Z) // stuck in a loop
+        break;
+
+      lc = INCREASED_Z;
+      transl.z() += 1;
+    }
+    else if(canonical_p.x() >= domain.xmax())
+    {
+      if(lc == INCREASED_X) // stuck in a loop
+        break;
+
+      lc = DECREASED_X;
+      transl.x() -= 1;
+    }
+    else if(canonical_p.y() >= domain.ymax())
+    {
+      if(lc == INCREASED_Y) // stuck in a loop
+        break;
+
+      lc = DECREASED_Y;
+      transl.y() -= 1;
+    }
+    else if(canonical_p.z() >= domain.zmax())
+    {
+      if(lc == INCREASED_Z) // stuck in a loop
+        break;
+
+      lc = DECREASED_Z;
+      transl.z() -= 1;
+    }
+    else
+    {
+      break;
+    }
+  }
+
+  if(canonical_p.x() < domain.xmin() || canonical_p.x() >= domain.xmax() ||
+     canonical_p.y() < domain.ymin() || canonical_p.y() >= domain.ymax() ||
+     canonical_p.z() < domain.zmin() || canonical_p.z() >= domain.zmax())
+  {
+    canonical_p = clamp_to_canonical_domain(canonical_p, gt);
   }
 
   CGAL_postcondition( !(canonical_p.x() < domain.xmin()) && (canonical_p.x() < domain.xmax()));
@@ -234,7 +253,7 @@ robust_canonicalize_point(const typename Gt_::Point_3& p,
 /// instance of the same weighted point that lives inside the base domain
 template <typename Gt_>
 typename Gt_::Weighted_point_3
-robust_canonicalize_point(const typename Gt_::Weighted_point_3& wp,
+construct_canonical_point(const typename Gt_::Weighted_point_3& wp,
                           const Gt_& gt)
 {
   typedef Gt_                                         Geom_traits;
@@ -252,7 +271,7 @@ robust_canonicalize_point(const typename Gt_::Weighted_point_3& wp,
   typename Geom_traits::Construct_weighted_point_3 cwp = gt.construct_weighted_point_3_object();
 
   const Bare_point& bp = cp(wp);
-  Bare_point canonical_point = robust_canonicalize_point(bp, gt);
+  Bare_point canonical_point = construct_canonical_point(bp, gt);
 
   return cwp(canonical_point, cw(wp));
 }