diff --git a/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp b/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp index 2d86e6a59d0..8e72bc854a7 100644 --- a/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp +++ b/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 } diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_function_wrapper.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_function_wrapper.h index bfb4b4e8716..040bfd7fff9 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_function_wrapper.h +++ b/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: diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 6d8201c773d..387e4125566 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/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 // = 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 pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits()); - std::pair pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits()); + std::pair pp_p = P3T3::internal::construct_periodic_point(p, geom_traits()); + std::pair 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 pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits()); + std::pair 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 pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits()); - std::pair pp_q = P3T3::internal::construct_periodic_point(cp(q), used_exact, geom_traits()); - std::pair pp_r = P3T3::internal::construct_periodic_point(cp(r), used_exact, geom_traits()); + std::pair pp_p = P3T3::internal::construct_periodic_point(p, geom_traits()); + std::pair pp_q = P3T3::internal::construct_periodic_point(cp(q), geom_traits()); + std::pair 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, diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h index e950cd7e4cf..957b73bec46 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h +++ b/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. diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h index 516c3e7434d..6c9c319c71e 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h +++ b/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()); } // --------------------------------------------------------------------------- diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h index 78e304e2949..e6b1ac20759 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h +++ b/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 std::pair 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,109 +66,46 @@ 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 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_::Point_3 -constrain_to_canonical_domain(const typename Gt_::Point_3& p, - const Gt_& gt) +clamp_to_canonical_domain(const typename Gt_::Point_3& p, + const Gt_& gt) { typedef Gt_ Geom_traits; typedef typename Geom_traits::FT FT; @@ -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_::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; + typedef Gt_ Geom_traits; + 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 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 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_::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)); }