From f509cbe8312e4d5bb5fe17e7b6db26d5b3a7a34d Mon Sep 17 00:00:00 2001 From: Aymeric PELLE Date: Thu, 26 Jun 2014 16:42:15 +0200 Subject: [PATCH] Fix the problem of holes on surface caused by bad analysis of dual facets. Add redefined predicats in Periodic_labeled_mesh_domain_3. (Do_intersect_surface and Construct_intersection.) Mesh_3 has new protected member functions allowing its heirs to access the labeling function and the error bound. Some private member functions required in Periodic_labeled_mesh_domain_3 are set protected. (Should be moved in a small feature modifying only Mesh_3?) --- Mesh_3/include/CGAL/Labeled_mesh_domain_3.h | 4 +- .../CGAL/Periodic_labeled_mesh_domain_3.h | 286 ++++++++++++++++++ 2 files changed, 289 insertions(+), 1 deletion(-) diff --git a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h index 8b051750d6e..912265f64c3 100644 --- a/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Labeled_mesh_domain_3.h @@ -422,7 +422,7 @@ public: // ----------------------------------- -private: +protected: /// Returns Surface_patch_index from \c i and \c j Surface_patch_index make_surface_index(const Subdomain_index i, const Subdomain_index j) const @@ -466,6 +466,8 @@ private: protected: /// Returns bounding box const Iso_cuboid_3& bounding_box() const { return bbox_; } + const Function& labeling_function() const { return function_; } + FT squared_error_bound_value() const { return squared_error_bound_; } private: /// The function which answers subdomain queries diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_labeled_mesh_domain_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_labeled_mesh_domain_3.h index eda3b0c47b3..861cf96ee17 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_labeled_mesh_domain_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_labeled_mesh_domain_3.h @@ -29,6 +29,15 @@ public: typedef typename Base::Sphere_3 Sphere_3; typedef typename Base::Bbox_3 Bbox_3; typedef typename Base::Iso_cuboid_3 Iso_cuboid_3; + typedef typename Base::Surface_patch Surface_patch; + typedef typename Base::Surface_patch_index Surface_patch_index; + typedef typename Base::Index Index; + typedef typename Base::Intersection Intersection; + typedef typename Base::Subdomain_index Subdomain_index; + typedef typename Base::Segment_3 Segment_3; + typedef typename Base::Ray_3 Ray_3; + typedef typename Base::Line_3 Line_3; + typedef typename Base::Point_3 Point_3; Periodic_labeled_mesh_domain_3(const Function& f, const Iso_cuboid_3& bbox, @@ -39,6 +48,283 @@ public: const Iso_cuboid_3& periodic_bounding_box() const { return Base::bounding_box(); } + struct Do_intersect_surface + { + Do_intersect_surface(const Periodic_labeled_mesh_domain_3& domain) + : r_domain_(domain) {} + + Surface_patch operator()(const Segment_3& s) const + { + return this->operator()(s.source(), s.target()); + } + + Surface_patch operator()(const Ray_3& r) const + { + return clip_to_segment(r); + } + + Surface_patch operator()(const Line_3& l) const + { + return clip_to_segment(l); + } + + private: + /// Returns true if points \c a & \c b do not belong to the same subdomain + /// \c index is set to the surface index of subdomains f(a), f(b) + Surface_patch operator()(const Point_3& a, const Point_3& b) const + { + // If f(a) != f(b), then [a,b] intersects some surface. Here we consider + // [a,b] intersects surface_patch labelled (or ). + // It may be false, further rafinement will improve precision + + Iso_cuboid_3 pbb = r_domain_.periodic_bounding_box(); + FT dimension [3] = { pbb.xmax()-pbb.xmin(), pbb.ymax()-pbb.ymin(), pbb.zmax()-pbb.zmin() }; + FT a_t [3] = { a.x(), a.y(), a.z() }; + FT b_t [3] = { b.x(), b.y(), b.z() }; + a_t[0] -= pbb.xmin(); + a_t[1] -= pbb.ymin(); + a_t[2] -= pbb.zmin(); + b_t[0] -= pbb.xmin(); + b_t[1] -= pbb.ymin(); + b_t[2] -= pbb.zmin(); + int o1 [3] = { a_t[0] / dimension[0], a_t[1] / dimension[1], a_t[2] / dimension[2] }; + int o2 [3] = { a_t[0] / dimension[0], a_t[1] / dimension[1], a_t[2] / dimension[2] }; + + unsigned i = 0; + while (i < 3) + { + if (abs(o1[i] - o2[i]) == 1) + break; + ++i; + } + if (i == 3) + { + Subdomain_index value_a = r_domain_.labeling_function()(a); + Subdomain_index value_b = r_domain_.labeling_function()(b); + if ( value_a != value_b ) + return Surface_patch(r_domain_.make_surface_index(value_a, value_b)); + return Surface_patch(); + } + + FT min_o = static_cast(std::min(o1[i], o2[i])); + FT a_min [3] = { a.x(), a.y(), a.z() }; + FT b_min [3] = { b.x(), b.y(), b.z() }; + a_min[i] -= (dimension[i] * min_o); + b_min[i] -= (dimension[i] * min_o); + + Subdomain_index value_a = r_domain_.labeling_function()(Point_3(a_min[0], a_min[1], a_min[2])); + Subdomain_index value_b = r_domain_.labeling_function()(Point_3(b_min[0], b_min[1], b_min[2])); + if ( value_a != value_b ) + return Surface_patch(r_domain_.make_surface_index(value_a, value_b)); + + FT max_o = static_cast(std::max(o1[i], o2[i])); + FT a_max [3] = { a.x(), a.y(), a.z() }; + FT b_max [3] = { b.x(), b.y(), b.z() }; + a_max[i] -= (dimension[i] * max_o); + b_max[i] -= (dimension[i] * max_o); + + value_a = r_domain_.labeling_function()(Point_3(a_max[0], a_max[1], a_max[2])); + value_b = r_domain_.labeling_function()(Point_3(b_max[0], b_max[1], b_max[2])); + if ( value_a != value_b ) + { + std::cout << "Fix works!" << std::endl; + return Surface_patch(r_domain_.make_surface_index(value_a, value_b)); + } + + return Surface_patch(); + } + + /** + * Clips \c query to a segment \c s, and call operator()(s) + */ + template + Surface_patch clip_to_segment(const Query& query) const + { + typename cpp11::result_of::type + clipped = CGAL::intersection(query, r_domain_.bounding_box()); + + if(clipped) +#if CGAL_INTERSECTION_VERSION > 1 + if(const Segment_3* s = boost::get(&*clipped)) + return this->operator()(*s); +#else + if(const Segment_3* s = object_cast(&clipped)) + return this->operator()(*s); +#endif + + return Surface_patch(); + } + + private: + const Periodic_labeled_mesh_domain_3& r_domain_; + }; + + /// Returns Do_intersect_surface object + Do_intersect_surface do_intersect_surface_object() const + { + return Do_intersect_surface(*this); + } + + struct Construct_intersection + { + Construct_intersection(const Periodic_labeled_mesh_domain_3& domain) + : r_domain_(domain) {} + + Intersection operator()(const Segment_3& s) const + { +#ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3 + CGAL_precondition(r_domain_.do_intersect_surface_object()(s)); +#endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3 + return this->operator()(s.source(),s.target()); + } + + Intersection operator()(const Ray_3& r) const + { + return clip_to_segment(r); + } + + Intersection operator()(const Line_3& l) const + { + return clip_to_segment(l); + } + + private: + /** + * Returns a point in the intersection of [a,b] with the surface + * \c a must be the source point, and \c b the out point. It's important + * because it drives bisection cuts. + * Indeed, the returned point is the first intersection from \c [a,b] + * with a subdomain surface. + */ + Intersection operator()(const Point_3& a, const Point_3& b) const + { + // Functors + typename BGT::Compute_squared_distance_3 squared_distance = + BGT().compute_squared_distance_3_object(); + typename BGT::Construct_midpoint_3 midpoint = + BGT().construct_midpoint_3_object(); + + Iso_cuboid_3 pbb = r_domain_.periodic_bounding_box(); + FT dimension [3] = { pbb.xmax()-pbb.xmin(), pbb.ymax()-pbb.ymin(), pbb.zmax()-pbb.zmin() }; + int o1 [3] = { a.x() / dimension[0], a.y() / dimension[1], a.z() / dimension[2] }; + int o2 [3] = { b.x() / dimension[0], b.y() / dimension[1], b.z() / dimension[2] }; + + unsigned i = 0; + while (i < 3) + { + if (abs(o1[i] - o2[i]) == 1) + break; + ++i; + } + if (i == 3) + i = 0; + + FT min_o = static_cast(std::min(o1[i], o2[i])); + FT a_min [3] = { a.x(), a.y(), a.z() }; + FT b_min [3] = { b.x(), b.y(), b.z() }; + a_min[i] -= (dimension[i] * min_o); + b_min[i] -= (dimension[i] * min_o); + + // Non const points + Point_3 p1(a_min[0], a_min[1], a_min[2]); + Point_3 p2(b_min[0], b_min[1], b_min[2]); + Point_3 mid = midpoint(p1, p2); + + // Cannot be const: those values are modified below. + Subdomain_index value_at_p1 = r_domain_.labeling_function()(p1); + Subdomain_index value_at_p2 = r_domain_.labeling_function()(p2); + Subdomain_index value_at_mid = r_domain_.labeling_function()(mid,true); + + // If both extremities are in the same subdomain, + // there is no intersection. + // This should not happen... + if( value_at_p1 == value_at_p2 ) + { + FT max_o = static_cast(std::max(o1[i], o2[i])); + FT a_max [3] = { a.x(), a.y(), a.z() }; + FT b_max [3] = { b.x(), b.y(), b.z() }; + a_max[i] -= (dimension[i] * max_o); + b_max[i] -= (dimension[i] * max_o); + p1 = Point_3(a_max[0], a_max[1], a_max[2]); + p2 = Point_3(b_max[0], b_max[1], b_max[2]); + mid = midpoint(p1, p2); + + value_at_p1 = r_domain_.labeling_function()(p1); + value_at_p2 = r_domain_.labeling_function()(p2); + value_at_mid = r_domain_.labeling_function()(mid,true); + + if( value_at_p1 == value_at_p2 ) + return Intersection(); + } + + // Construct the surface patch index and index from the values at 'a' + // and 'b'. Even if the bissection find out a different pair of + // values, the reported index will be constructed from the initial + // values. + const Surface_patch_index sp_index = + r_domain_.make_surface_index(value_at_p1, value_at_p2); + const Index index = r_domain_.index_from_surface_patch_index(sp_index); + + // Else lets find a point (by bisection) + // Bisection ends when the point is near than error bound from surface + while(true) + { + // If the two points are enough close, then we return midpoint + if ( squared_distance(p1, p2) < r_domain_.squared_error_bound_value() ) + { + CGAL_assertion(value_at_p1 != value_at_p2); + return Intersection(mid, index, 2); + } + + // Else we must go on + // Here we consider that p1(a) is the source point. Thus, we keep p1 and + // change p2 if f(p1)!=f(p2). + // That allows us to find the first intersection from a of [a,b] with + // a surface. + if ( value_at_p1 != value_at_mid ) + { + p2 = mid; + value_at_p2 = value_at_mid; + } + else + { + p1 = mid; + value_at_p1 = value_at_mid; + } + + mid = midpoint(p1, p2); + value_at_mid = r_domain_.labeling_function()(mid,true); + } + } + + /// Clips \c query to a segment \c s, and call operator()(s) + template + Intersection clip_to_segment(const Query& query) const + { + typename cpp11::result_of::type + clipped = CGAL::intersection(query, r_domain_.bounding_box()); + + if(clipped) +#if CGAL_INTERSECTION_VERSION > 1 + if(const Segment_3* s = boost::get(&*clipped)) + return this->operator()(*s); +#else + if(const Segment_3* s = object_cast(&clipped)) + return this->operator()(*s); +#endif + + return Intersection(); + } + + private: + const Periodic_labeled_mesh_domain_3& r_domain_; + }; + + Construct_intersection construct_intersection_object() const + { + return Construct_intersection(*this); + } + private: // Disabled copy constructor & assignment operator Periodic_labeled_mesh_domain_3(const Periodic_labeled_mesh_domain_3&);