From 19e5e4a16e1331ade755b4c29b2b2fe8a2b4264c Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 2 Mar 2021 12:25:28 +0000 Subject: [PATCH 01/27] Use flaat_set and put it outside the loop --- .../Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h index 51b6ebba6830..9e704870c183 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -177,14 +178,16 @@ bool is_polygon_soup_a_polygon_mesh(const PolygonRange& polygons) //check there is no duplicated ordered edge, and //check there is no polygon with twice the same vertex std::set > edge_set; + boost::container::flat_set polygon_vertices; V_ID max_id = 0; + for(const Polygon& polygon : polygons) { std::size_t nb_edges = boost::size(polygon); if(nb_edges < 3) return false; - std::set polygon_vertices; + polygon_vertices.clear(); V_ID prev = *std::prev(boost::end(polygon)); for(V_ID id : polygon) { From 10f454d98384dfacb2cd5a06b1f7df6a7ff6ea4a Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 2 Mar 2021 13:09:45 +0000 Subject: [PATCH 02/27] reserve(8) for incident faces to a vertex --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index e31e91db8cdb..7f438d4d5bfa 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -506,6 +506,7 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits incident_faces; + incident_faces.reserve(8); for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh)) { if(f == boost::graph_traits::null_face()) From 5e2a580551527336ff8dd13357b1465ad4334d23 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 2 Mar 2021 13:17:53 +0000 Subject: [PATCH 03/27] map -> unordered_map --- .../include/CGAL/Polygon_mesh_processing/compute_normal.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index 7f438d4d5bfa..bfb8f831238e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -30,6 +30,7 @@ #include #include #include +#include #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP # ifndef CGAL_PMP_COMPUTE_NORMAL_DEBUG @@ -670,7 +671,7 @@ compute_vertex_normal(typename boost::graph_traits::vertex_descript VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pmesh)); - typedef std::map Face_vector_map; + typedef std::unordered_map Face_vector_map; typedef boost::associative_property_map Default_map; typedef typename internal_np::Lookup_named_param_def Date: Tue, 2 Mar 2021 14:40:22 +0000 Subject: [PATCH 04/27] Important gain for Euler::add_face() --- BGL/include/CGAL/boost/graph/Euler_operations.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/BGL/include/CGAL/boost/graph/Euler_operations.h b/BGL/include/CGAL/boost/graph/Euler_operations.h index 8d5268da6e0d..beed0ecfb39d 100644 --- a/BGL/include/CGAL/boost/graph/Euler_operations.h +++ b/BGL/include/CGAL/boost/graph/Euler_operations.h @@ -737,9 +737,9 @@ add_face(const VertexRange& vr, Graph& g) patch_start, patch_end; // cache for set_next and vertex' set_halfedge typedef std::pair NextCacheEntry; - typedef std::vector NextCache; + typedef boost::container::small_vector NextCache; NextCache next_cache; - next_cache.reserve(3 * n); + //next_cache.reserve(3 * n); Unfortunately small_vector has no reserve // re-link patches if necessary for (unsigned int i = 0, ii = 1; i Date: Tue, 2 Mar 2021 21:30:36 +0000 Subject: [PATCH 05/27] small_vector has reserve() --- BGL/include/CGAL/boost/graph/Euler_operations.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BGL/include/CGAL/boost/graph/Euler_operations.h b/BGL/include/CGAL/boost/graph/Euler_operations.h index beed0ecfb39d..d26d6a4bcebe 100644 --- a/BGL/include/CGAL/boost/graph/Euler_operations.h +++ b/BGL/include/CGAL/boost/graph/Euler_operations.h @@ -739,7 +739,7 @@ add_face(const VertexRange& vr, Graph& g) typedef std::pair NextCacheEntry; typedef boost::container::small_vector NextCache; NextCache next_cache; - //next_cache.reserve(3 * n); Unfortunately small_vector has no reserve + next_cache.reserve(3 * n); // re-link patches if necessary for (unsigned int i = 0, ii = 1; i Date: Tue, 2 Mar 2021 21:56:54 +0000 Subject: [PATCH 06/27] Replace geometric by combinatorial test --- .../internal/Isotropic_remeshing/remesh_impl.h | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 1a0fa1d3b8f6..ce1636b62d08 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -1393,6 +1393,7 @@ namespace internal { bool collapse_would_invert_face(const halfedge_descriptor& h) const { + vertex_descriptor tv = target(h, mesh_); typename boost::property_traits::reference s = get(vpmap_, source(h, mesh_)); //s for source typename boost::property_traits::reference @@ -1409,16 +1410,21 @@ namespace internal { if (face(hd, mesh_) == boost::graph_traits::null_face()) continue; + vertex_descriptor tnhd = target(next(hd, mesh_), mesh_); + vertex_descriptor tnnhd = target(next(next(hd, mesh_), mesh_), mesh_); typename boost::property_traits::reference - p = get(vpmap_, target(next(hd, mesh_), mesh_)); + p = get(vpmap_, tnhd); typename boost::property_traits::reference - q = get(vpmap_, target(next(next(hd, mesh_), mesh_), mesh_)); + q = get(vpmap_, tnnhd); #ifdef CGAL_PMP_REMESHING_DEBUG CGAL_assertion((Triangle_3(t, p, q).is_degenerate()) == GeomTraits().collinear_3_object()(t, p, q)); #endif + if((tv == tnnhd) || (tv == tnhd)) + continue; + if ( GeomTraits().collinear_3_object()(s, p, q) || GeomTraits().collinear_3_object()(t, p, q)) continue; From 2ac977b7c5df414979a515033bb07c0d16cc02de Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 13:21:24 +0000 Subject: [PATCH 07/27] Early exit in do_intersect of Sphere/Bbox_3 --- .../internal/Bbox_3_Sphere_3_do_intersect.h | 34 +++++++++++++++---- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h index b15270e52942..49eaf967ad5b 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h @@ -35,39 +35,59 @@ namespace internal { typedef typename K::Point_3 Point; FT d = FT(0); FT distance = FT(0); + FT sr = sphere.squared_radius(); + Point center = sphere.center(); if(center.x() < (FT)bbox.xmin()) { d = (FT)bbox.xmin() - center.x(); - distance += d * d; + d = square(d); + if(d > sr){ + return false; + } + distance = d; } else if(center.x() > (FT)bbox.xmax()) { d = center.x() - (FT)bbox.xmax(); - distance += d * d; + d = square(d); + if(d > sr){ + return false; + } + distance = d; } if(center.y() < (FT)bbox.ymin()) { d = (FT)bbox.ymin() - center.y(); - distance += d * d; + d = square(d); + if(d > sr){ + return false; + } + distance += d ; } else if(center.y() > (FT)bbox.ymax()) { d = center.y() - (FT)bbox.ymax(); - distance += d * d; + d = square(d); + if(d > sr){ + return false; + } + distance += d; } if(center.z() < (FT)bbox.zmin()) { d = (FT)bbox.zmin() - center.z(); - distance += d * d; + d = square(d); + distance += d; } else if(center.z() > (FT)bbox.zmax()) { d = center.z() - (FT)bbox.zmax(); - distance += d * d; + d = square(d); + distance += d; } // For unknown reason this causes a syntax error on VC2005 @@ -87,7 +107,7 @@ namespace internal { // } //} - return distance <= sphere.squared_radius(); + return distance <= sr; } template From 82e5b3a02d51ec2ce672bc00736799a80ba3c809 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 14:23:56 +0000 Subject: [PATCH 08/27] Do not collect in a vector per patch, but do the tests directly (@janetournois please double-check this commit) --- .../Isotropic_remeshing/remesh_impl.h | 39 ++++++------------- 1 file changed, 12 insertions(+), 27 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index ce1636b62d08..edb6f3e0dbc4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -43,6 +43,7 @@ #include #include #include +#include #include #include @@ -1885,8 +1886,8 @@ namespace internal { bool check_normals(const HalfedgeRange& hedges) const { std::size_t nb_patches = patch_id_to_index_map.size(); + std::vector > normal_per_patch(nb_patches,boost::none); - std::vector< std::vector > normals_per_patch(nb_patches); for(halfedge_descriptor hd : hedges) { Halfedge_status s = status(hd); @@ -1898,18 +1899,15 @@ namespace internal { if (n == CGAL::NULL_VECTOR) //for degenerate faces continue; Patch_id pid = get_patch_id(face(hd, mesh_)); - normals_per_patch[patch_id_to_index_map.at(pid)].push_back(n); - } - - //on each surface patch, - //check all normals have same orientation - for (std::size_t i=0; i < nb_patches; ++i) - { - const std::vector& normals = normals_per_patch[i]; - if (normals.empty()) continue; - - if (!check_orientation(normals)) - return false; + std::size_t index = patch_id_to_index_map.at(pid); + if(normal_per_patch[index]){ + const Vector_3& vec = *normal_per_patch[index]; + double dot = to_double(n * vec); + if (dot <= 0.){ + return false; + } + } + normal_per_patch[index] = boost::make_optional(n); } return true; } @@ -1922,20 +1920,7 @@ namespace internal { Vector_3 no = compute_normal(face(opposite(h, mesh_), mesh_)); return n * no > 0.; } - - bool check_orientation(const std::vector& normals) const - { - if (normals.size() < 2) - return true; - for (std::size_t i = 1; i < normals.size(); ++i)/*start at 1 on purpose*/ - { - double dot = to_double(normals[i - 1] * normals[i]); - if (dot <= 0.) - return false; - } - return true; - } - + public: const Triangle_list& input_triangles() const { return input_triangles_; From dad0287b69b5b347b1c687f8fc363b343bc42f73 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 14:34:56 +0000 Subject: [PATCH 09/27] Remove trailing whitespace --- .../internal/Isotropic_remeshing/remesh_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index edb6f3e0dbc4..d0bd9bf2a008 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -1920,7 +1920,7 @@ namespace internal { Vector_3 no = compute_normal(face(opposite(h, mesh_), mesh_)); return n * no > 0.; } - + public: const Triangle_list& input_triangles() const { return input_triangles_; From e7f18dff56952c7a3e4999eac4ae0fd14fff3eaa Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 15:00:59 +0000 Subject: [PATCH 10/27] No need for calling abs for a squared distance --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 92cbd3ebf3ff..88875d6d3117 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -245,18 +245,17 @@ class Do_intersect_3 } double double_tmp_result = (distance - ssr); - double max2 = CGAL::abs(ssr); - if ((max1 < 3.33558365626356687717e-147) || (max2 < 1.11261183279326254436e-293)){ + if ((max1 < 3.33558365626356687717e-147) || (ssr < 1.11261183279326254436e-293)){ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } - if ((max1 > 1.67597599124282407923e+153) || (max2 > 2.80889552322236673473e+306)){ + if ((max1 > 1.67597599124282407923e+153) || (ssr > 2.80889552322236673473e+306)){ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } - double eps = 1.99986535548615598560e-15 * (std::max) (max2, (max1 * max1)); + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); if (double_tmp_result > eps) return false; @@ -265,6 +264,7 @@ class Do_intersect_3 if (double_tmp_result < -eps) return true; else { + return true; // AF: just to test conservative approach CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } From 431c5b81322801ad6740c68987d65d2bca1a7e3a Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 15:02:13 +0000 Subject: [PATCH 11/27] Comment something experimental that sneaked in --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 88875d6d3117..9c09fe2ecb33 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -264,7 +264,7 @@ class Do_intersect_3 if (double_tmp_result < -eps) return true; else { - return true; // AF: just to test conservative approach + // return true; // AF: just to test conservative approach CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } From 7785a9c2b142e2034811b7bd8a47ae2864581ef0 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 15:25:30 +0000 Subject: [PATCH 12/27] Use certainly() --- .../internal/Bbox_3_Sphere_3_do_intersect.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h index 49eaf967ad5b..7f74ddf1c8ef 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h @@ -43,7 +43,7 @@ namespace internal { { d = (FT)bbox.xmin() - center.x(); d = square(d); - if(d > sr){ + if(certainly(d > sr)){ return false; } distance = d; @@ -52,7 +52,7 @@ namespace internal { { d = center.x() - (FT)bbox.xmax(); d = square(d); - if(d > sr){ + if(certainly(d > sr)){ return false; } distance = d; @@ -62,7 +62,7 @@ namespace internal { { d = (FT)bbox.ymin() - center.y(); d = square(d); - if(d > sr){ + if(certainly(d > sr)){ return false; } distance += d ; @@ -71,7 +71,7 @@ namespace internal { { d = center.y() - (FT)bbox.ymax(); d = square(d); - if(d > sr){ + if(certainly(d > sr)){ return false; } distance += d; From 4fab8430dad1506bc6c7eb36c0f4dbb993a81e55 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 3 Mar 2021 21:14:39 +0000 Subject: [PATCH 13/27] Add static filter for Equal_3::operator()(Vector_3, Null_vector) --- .../CGAL/internal/Static_filters/Equal_3.h | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Equal_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Equal_3.h index 102199d56335..d4e3c3099ed1 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Equal_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Equal_3.h @@ -82,6 +82,24 @@ class Equal_3 return Base::operator()(p, q); } + + result_type operator()(const Vector_3 &p, const Null_vector &q) const + { + CGAL_BRANCH_PROFILER(std::string("semi-static attempts/calls to : ") + + std::string(CGAL_PRETTY_FUNCTION), tmp); + + Get_approx get_approx; // Identity functor for all points + // but lazy points + double px, py, pz; + + if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) && + fit_in_double(get_approx(p).z(), pz) ) + { + CGAL_BRANCH_PROFILER_BRANCH(tmp); + return px == 0 && py == 0 && pz == 0; + } + return Base::operator()(p, q); + } }; // end class Equal_3 } // end namespace Static_filters_predicates From 2e3bfa8743a9721abd4be13f3e271af5f2f0eb7c Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 4 Mar 2021 09:27:05 +0000 Subject: [PATCH 14/27] Reduce calls to target @sloriot please double check the correctness --- .../CGAL/Polygon_mesh_processing/compute_normal.h | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h index bfb8f831238e..0f5840f6667a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h @@ -91,12 +91,16 @@ void sum_normals(const PM& pmesh, halfedge_descriptor he = halfedge(f, pmesh); vertex_descriptor v = source(he, pmesh); + vertex_descriptor the = target(he,pmesh); + he = next(he, pmesh); + vertex_descriptor tnhe = target(he,pmesh); + const Point_ref pv = get(vpmap, v); - while(v != target(next(he, pmesh), pmesh)) + while(v != tnhe) { - const Point_ref pvn = get(vpmap, target(he, pmesh)); - const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); + const Point_ref pvn = get(vpmap, the); + const Point_ref pvnn = get(vpmap, tnhe); const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); @@ -107,7 +111,9 @@ void sum_normals(const PM& pmesh, sum = traits.construct_sum_of_vectors_3_object()(sum, n); + the = tnhe; he = next(he, pmesh); + tnhe = target(he,pmesh); } } From 45265a7e342955f498a062ffb9845c532efd119b Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 4 Mar 2021 10:46:46 +0000 Subject: [PATCH 15/27] Replace vector with optional by two vectors. No idea what is better yet --- .../internal/Isotropic_remeshing/remesh_impl.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index d0bd9bf2a008..e2399b115e72 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -1886,7 +1886,9 @@ namespace internal { bool check_normals(const HalfedgeRange& hedges) const { std::size_t nb_patches = patch_id_to_index_map.size(); - std::vector > normal_per_patch(nb_patches,boost::none); + //std::vector > normal_per_patch(nb_patches,boost::none); + std::vector initialized(nb_patches,false); + std::vector normal_per_patch(nb_patches); for(halfedge_descriptor hd : hedges) { @@ -1900,14 +1902,17 @@ namespace internal { continue; Patch_id pid = get_patch_id(face(hd, mesh_)); std::size_t index = patch_id_to_index_map.at(pid); - if(normal_per_patch[index]){ - const Vector_3& vec = *normal_per_patch[index]; + //if(normal_per_patch[index]){ + if(initialized[index]){ + const Vector_3& vec = normal_per_patch[index]; double dot = to_double(n * vec); if (dot <= 0.){ return false; } } - normal_per_patch[index] = boost::make_optional(n); + //normal_per_patch[index] = boost::make_optional(n); + normal_per_patch[index] = n; + initialized[index] = true; } return true; } From b1c2dd8db67f5149c6080a91359ad3a81d38b5bc Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Sun, 7 Mar 2021 18:25:18 +0000 Subject: [PATCH 16/27] WIP: early exit in the static filter --- AABB_tree/include/CGAL/AABB_traits.h | 20 +++ .../internal/Static_filters/Do_intersect_3.h | 132 +++++++++++++++++- .../isotropic_remeshing_example.cpp | 23 ++- 3 files changed, 169 insertions(+), 6 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 613d08c174bf..1a3fdde10551 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -423,6 +423,26 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } + /* + template <> + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& pr, const Point& bound) const + { + return GeomTraits().do_intersect_3_object() + (GeomTraits().construct_sphere_3_object() + (p, GeomTraits().compute_squared_distance_3_object()(p, bound)).bbox(), pr)? + CGAL::SMALLER : CGAL::LARGER; + } + + template <> + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& pr, const FT& sq_distance) const + { + return GeomTraits().do_intersect_3_object() + (GeomTraits().construct_sphere_3_object()(p, sq_distance).bbox(), + pr) ? + CGAL::SMALLER : + CGAL::LARGER; + } + */ }; Closest_point closest_point_object() const {return Closest_point(*this);} diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 9c09fe2ecb33..8c899a32c181 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -192,6 +192,10 @@ class Do_intersect_3 { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); + if ((ssr < 1.11261183279326254436e-293) || (ssr > 2.80889552322236673473e+306)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } double distance = 0; double max1 = 0; @@ -201,6 +205,25 @@ class Do_intersect_3 max1 = bxmin_scx; distance = square(bxmin_scx); +#ifdef EARLY + { + double double_tmp_result = (distance - ssr); + + if (max1 < 3.33558365626356687717e-147){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if (max1 > 1.67597599124282407923e+153){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } else if(scx > bxmax) { @@ -208,8 +231,31 @@ class Do_intersect_3 max1 = scx_bxmax; distance = square(scx_bxmax); +#ifdef EARLY + { + double double_tmp_result = (distance - ssr); + + if (max1 < 3.33558365626356687717e-147){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if (max1 > 1.67597599124282407923e+153){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } + + + + if(scy < bymin) { double bymin_scy = bymin - scy; @@ -217,6 +263,25 @@ class Do_intersect_3 max1 = bymin_scy; distance += square(bymin_scy); +#ifdef EARLY + { + double double_tmp_result = (distance - ssr); + + if (max1 < 3.33558365626356687717e-147){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if ((max1 > 1.67597599124282407923e+153)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } else if(scy > bymax) { @@ -225,6 +290,26 @@ class Do_intersect_3 max1 = scy_bymax; distance += square(scy_bymax); +#ifdef EARLY + + { + double double_tmp_result = (distance - ssr); + + if ((max1 < 3.33558365626356687717e-147)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if ((max1 > 1.67597599124282407923e+153)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } if(scz < bzmin) @@ -234,6 +319,25 @@ class Do_intersect_3 max1 = bzmin_scz; distance += square(bzmin_scz); +#ifdef EARLY + { + double double_tmp_result = (distance - ssr); + + if ((max1 < 3.33558365626356687717e-147)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if ((max1 > 1.67597599124282407923e+153)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } else if(scz > bzmax) { @@ -242,19 +346,39 @@ class Do_intersect_3 max1 = scz_bzmax; distance += square(scz_bzmax); +#ifdef EARLY + { + double double_tmp_result = (distance - ssr); + + if ((max1 < 3.33558365626356687717e-147)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + if ((max1 > 1.67597599124282407923e+153)){ + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); + } + + double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps) + return false; + } +#endif } double double_tmp_result = (distance - ssr); - if ((max1 < 3.33558365626356687717e-147) || (ssr < 1.11261183279326254436e-293)){ +#ifndef EARLY + if (max1 < 3.33558365626356687717e-147){ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } - if ((max1 > 1.67597599124282407923e+153) || (ssr > 2.80889552322236673473e+306)){ + if (max1 > 1.67597599124282407923e+153){ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } - +#endif double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); if (double_tmp_result > eps) @@ -264,7 +388,7 @@ class Do_intersect_3 if (double_tmp_result < -eps) return true; else { - // return true; // AF: just to test conservative approach + return true; // AF: just to test conservative approach CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp index c5811a87d2b2..bb246a9393c9 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp @@ -1,6 +1,7 @@ +#define EARLY #include #include - +#include #include #include #include @@ -9,6 +10,7 @@ #include #include +#include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Surface_mesh Mesh; @@ -35,6 +37,9 @@ int main(int argc, char* argv[]) { const char* filename = (argc > 1) ? argv[1] : "data/pig.off"; + CGAL::Timer t; + t.start(); + Mesh mesh; if(!PMP::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) { @@ -42,7 +47,10 @@ int main(int argc, char* argv[]) return 1; } - double target_edge_length = 0.04; + std::cout << "Read: " << t.time() << " sec" << std::endl; + t.reset(); + + double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 1.0; unsigned int nb_iter = 3; std::cout << "Split border..."; @@ -53,6 +61,9 @@ int main(int argc, char* argv[]) std::cout << "done." << std::endl; + std::cout << "took " << t.time() << " sec" << std::endl; + t.reset(); + std::cout << "Start remeshing of " << filename << " (" << num_faces(mesh) << " faces)..." << std::endl; @@ -62,5 +73,13 @@ int main(int argc, char* argv[]) std::cout << "Remeshing done." << std::endl; + std::cout << "took " << t.time() << " sec" << std::endl; + t.reset(); + + std::ofstream out("out.off"); + out.precision(17); + //out << mesh << std::endl; + out.close(); + return 0; } From 643810f310f4df42c48b68fe410560866f8db006 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 9 Mar 2021 13:15:49 +0000 Subject: [PATCH 17/27] WIP for not hard coding conservative. Todo: Use of Has_filtered --- AABB_tree/include/CGAL/AABB_traits.h | 23 +++++-------------- .../internal/Static_filters/Do_intersect_3.h | 8 +++++-- 2 files changed, 12 insertions(+), 19 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 1a3fdde10551..6bcf44a18afc 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -414,35 +414,24 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - template - CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const - { - return GeomTraits().do_intersect_3_object() - (GeomTraits().construct_sphere_3_object()(p, sq_distance), - pr) ? - CGAL::SMALLER : - CGAL::LARGER; - } - /* - template <> - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& pr, const Point& bound) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() - (p, GeomTraits().compute_squared_distance_3_object()(p, bound)).bbox(), pr)? + (p, GeomTraits().compute_squared_distance_3_object()(p, bound)), bb,true)? CGAL::SMALLER : CGAL::LARGER; } - template <> - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& pr, const FT& sq_distance) const + template + CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const { return GeomTraits().do_intersect_3_object() - (GeomTraits().construct_sphere_3_object()(p, sq_distance).bbox(), + (GeomTraits().construct_sphere_3_object()(p, sq_distance), pr) ? CGAL::SMALLER : CGAL::LARGER; } - */ + }; Closest_point closest_point_object() const {return Closest_point(*this);} diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index 8c899a32c181..d09972a871a7 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -168,9 +168,11 @@ class Do_intersect_3 return this->operator()(s, b); } + // The parameter overestimate is used to avoid a filter failure in AABB_tree::closest_point() result_type - operator()(const Sphere_3 &s, const Bbox_3& b) const + operator()(const Sphere_3 &s, const Bbox_3& b, bool overestimate = false) const { + std::cout << "." << std::flush; CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") + std::string(CGAL_PRETTY_FUNCTION), tmp); @@ -388,7 +390,9 @@ class Do_intersect_3 if (double_tmp_result < -eps) return true; else { - return true; // AF: just to test conservative approach + if(overestimate){ + return true; + } CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } From 9274503f125ad0082ca101e0f2055a9ce19b9276 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 9 Mar 2021 13:18:50 +0000 Subject: [PATCH 18/27] remove debug code --- .../include/CGAL/internal/Static_filters/Do_intersect_3.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index d09972a871a7..f5149d93f156 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -172,7 +172,6 @@ class Do_intersect_3 result_type operator()(const Sphere_3 &s, const Bbox_3& b, bool overestimate = false) const { - std::cout << "." << std::flush; CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") + std::string(CGAL_PRETTY_FUNCTION), tmp); From cea8cca6af0be8e1bfc8974b6540a78a008fe02f Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 9 Mar 2021 15:45:18 +0000 Subject: [PATCH 19/27] WIP: less computations of degree --- .../Isotropic_remeshing/remesh_impl.h | 50 +++++++++++++++++-- 1 file changed, 45 insertions(+), 5 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index e2399b115e72..551c1294cfaa 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -826,6 +827,22 @@ namespace internal { #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << "Equalize valences..." << std::endl; #endif + +#define VALE +#ifdef VALE + typedef typename boost::property_map >::type Vertex_degree; + Vertex_degree degree = get(CGAL::dynamic_vertex_property_t(), mesh_); + + for(vertex_descriptor v : vertices(mesh_)){ + put(degree,v,0); + } + for(halfedge_descriptor h : halfedges(mesh_)) + { + vertex_descriptor t = target(h, mesh_); + put(degree, t, get(degree,t)+1); + } +#endif + unsigned int nb_flips = 0; for(edge_descriptor e : edges(mesh_)) { @@ -838,11 +855,17 @@ namespace internal { vertex_descriptor vb = target(he, mesh_); vertex_descriptor vc = target(next(he, mesh_), mesh_); vertex_descriptor vd = target(next(opposite(he, mesh_), mesh_), mesh_); - - int vva = valence(va), tvva = target_valence(va); - int vvb = valence(vb), tvvb = target_valence(vb); - int vvc = valence(vc), tvvc = target_valence(vc); - int vvd = valence(vd), tvvd = target_valence(vd); +#ifdef VALE + int vva = get(degree,va), tvva = target_valence(va); + int vvb = get(degree, vb), tvvb = target_valence(vb); + int vvc = get(degree,vc), tvvc = target_valence(vc); + int vvd = get(degree,vd), tvvd = target_valence(vd); +#else + int vva = valence(va), tvva = target_valence(va); + int vvb = valence( vb), tvvb = target_valence(vb); + int vvc = valence(vc), tvvc = target_valence(vc); + int vvd = valence(vd), tvvd = target_valence(vd); +#endif int deviation_pre = CGAL::abs(vva - tvva) + CGAL::abs(vvb - tvvb) + CGAL::abs(vvc - tvvc) @@ -860,6 +883,12 @@ namespace internal { vvb -= 1; vvc += 1; vvd += 1; +#ifdef VALE + put(degree, va, vva); + put(degree, vb, vvb); + put(degree, vc, vvc); + put(degree, vd, vvd); +#endif ++nb_flips; #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS @@ -895,6 +924,17 @@ namespace internal { CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) ); CGAL_assertion( !get(ecmap_, edge(he, mesh_)) ); CGAL::Euler::flip_edge(he, mesh_); +#ifdef VALE + vva += 1; + vvb += 1; + vvc -= 1; + vvd -= 1; + + put(degree, va, vva); + put(degree, vb, vvb); + put(degree, vc, vvc); + put(degree, vd, vvd); +#endif --nb_flips; CGAL_assertion_code(Halfedge_status s3 = status(he)); From b5239161c46576daccc7d37b991c43974c656f27 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 10 Mar 2021 10:51:24 +0000 Subject: [PATCH 20/27] Check that the kernel has static filters --- AABB_tree/include/CGAL/AABB_traits.h | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 6bcf44a18afc..78959df857be 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -414,7 +414,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_true&) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() @@ -422,6 +422,19 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_false&) const + { + return GeomTraits().do_intersect_3_object() + (GeomTraits().construct_sphere_3_object() + (p, GeomTraits().compute_squared_distance_3_object()(p, bound)), bb)? + CGAL::SMALLER : CGAL::LARGER; + } + + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const + { + return (*this)(p, bb, bound, typename GeomTraits::Has_filtered_predicates_tag()); + } + template CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const { From ac708f6625e4a6448f88ac414cbf3c0fcc917c03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 10 Mar 2021 13:59:00 +0100 Subject: [PATCH 21/27] use nested class for dispatch --- AABB_tree/include/CGAL/AABB_traits.h | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 78959df857be..a022070a654a 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -24,6 +24,8 @@ #include #include #include +#include + #include @@ -414,7 +416,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_true&) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, const Tag_true&) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() @@ -422,7 +424,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_false&) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, const Tag_false&) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() @@ -430,10 +432,22 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const - { - return (*this)(p, bb, bound, typename GeomTraits::Has_filtered_predicates_tag()); - } + template ::value> + struct Use_conservative_static_filters + { + typedef CGAL::Tag_false type; + }; + + template + struct Use_conservative_static_filters + { + typedef CGAL::Boolean_tag type; + }; + + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const + { + return (*this)(p, bb, bound, typename Use_conservative_static_filters::type()); + } template CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const From 3614dee0fbe1fd4b414b32aef8c4b82f060d3756 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 10 Mar 2021 14:01:55 +0100 Subject: [PATCH 22/27] don't reinvent the wheel --- AABB_tree/include/CGAL/AABB_traits.h | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index a022070a654a..85c86355c570 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -432,21 +432,9 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - template ::value> - struct Use_conservative_static_filters - { - typedef CGAL::Tag_false type; - }; - - template - struct Use_conservative_static_filters - { - typedef CGAL::Boolean_tag type; - }; - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const { - return (*this)(p, bb, bound, typename Use_conservative_static_filters::type()); + return (*this)(p, bb, bound, CGAL::Boolean_tag::value>()); } template From 2b5cb2da9f08cc1fdbed23886c227dac16e5eaad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 10 Mar 2021 14:23:58 +0100 Subject: [PATCH 23/27] copy tags --- AABB_tree/include/CGAL/AABB_traits.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 85c86355c570..7cabbdae1df7 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -416,7 +416,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, const Tag_true&) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_true) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() @@ -424,7 +424,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } - CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, const Tag_false&) const + CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound, Tag_false) const { return GeomTraits().do_intersect_3_object() (GeomTraits().construct_sphere_3_object() @@ -434,7 +434,7 @@ class AABB_traits CGAL::Comparison_result operator()(const Point& p, const Bounding_box& bb, const Point& bound) const { - return (*this)(p, bb, bound, CGAL::Boolean_tag::value>()); + return (*this)(p, bb, bound, Boolean_tag::value>()); } template From e6c2e5c9499c4c94fb994e1de15f9daca51ce8af Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 10 Mar 2021 13:23:01 +0000 Subject: [PATCH 24/27] Simplify. todo: tighter bounds per case or factorization with a lambda --- .../internal/Static_filters/Do_intersect_3.h | 193 +++++++----------- 1 file changed, 77 insertions(+), 116 deletions(-) diff --git a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h index f5149d93f156..8f755a05965b 100644 --- a/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h @@ -178,18 +178,14 @@ class Do_intersect_3 Get_approx get_approx; // Identity functor for all points const Point_3& c = s.center(); - double scx, scy, scz, ssr, bxmin, bymin, bzmin, bxmax, bymax, bzmax; + double scx, scy, scz, ssr; + double bxmin = b.xmin() , bymin = b.ymin() , bzmin = b.zmin() , + bxmax = b.xmax() , bymax = b.ymax() , bzmax = b.zmax() ; if (fit_in_double(get_approx(c).x(), scx) && fit_in_double(get_approx(c).y(), scy) && fit_in_double(get_approx(c).z(), scz) && - fit_in_double(s.squared_radius(), ssr) && - fit_in_double(b.xmin(), bxmin) && - fit_in_double(b.ymin(), bymin) && - fit_in_double(b.zmin(), bzmin) && - fit_in_double(b.xmax(), bxmax) && - fit_in_double(b.ymax(), bymax) && - fit_in_double(b.zmax(), bzmax)) + fit_in_double(s.squared_radius(), ssr)) { CGAL_BRANCH_PROFILER_BRANCH_1(tmp); @@ -199,32 +195,30 @@ class Do_intersect_3 } double distance = 0; double max1 = 0; - + double double_tmp_result = 0; + double eps = 0; if(scx < bxmin) { double bxmin_scx = bxmin - scx; max1 = bxmin_scx; distance = square(bxmin_scx); -#ifdef EARLY - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if (max1 < 3.33558365626356687717e-147){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if (max1 > 1.67597599124282407923e+153){ + if( (max1 < 3.33558365626356687717e-147) || (max1 > 1.67597599124282407923e+153) ){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } + } - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; + if (double_tmp_result > eps){ + return false; } -#endif } else if(scx > bxmax) { @@ -232,169 +226,136 @@ class Do_intersect_3 max1 = scx_bxmax; distance = square(scx_bxmax); -#ifdef EARLY - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if (max1 < 3.33558365626356687717e-147){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if (max1 > 1.67597599124282407923e+153){ + if( (max1 < 3.33558365626356687717e-147) || (max1 > 1.67597599124282407923e+153)){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } + } - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; + if (double_tmp_result > eps){ + return false; } -#endif } - - - if(scy < bymin) { double bymin_scy = bymin - scy; - if(max1 < bymin_scy) + if(max1 < bymin_scy){ max1 = bymin_scy; + } distance += square(bymin_scy); -#ifdef EARLY - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if (max1 < 3.33558365626356687717e-147){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if ((max1 > 1.67597599124282407923e+153)){ + if( (max1 < 3.33558365626356687717e-147) || ((max1 > 1.67597599124282407923e+153)) ){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } + } - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; + if (double_tmp_result > eps){ + return false; } -#endif } else if(scy > bymax) { double scy_bymax = scy - bymax; - if(max1 < scy_bymax) + if(max1 < scy_bymax){ max1 = scy_bymax; - + } distance += square(scy_bymax); -#ifdef EARLY - - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if ((max1 < 3.33558365626356687717e-147)){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if ((max1 > 1.67597599124282407923e+153)){ + if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153)) ){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } + } - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; + if (double_tmp_result > eps){ + return false; } -#endif } + if(scz < bzmin) { double bzmin_scz = bzmin - scz; - if(max1 < bzmin_scz) + if(max1 < bzmin_scz){ max1 = bzmin_scz; - + } distance += square(bzmin_scz); -#ifdef EARLY - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if ((max1 < 3.33558365626356687717e-147)){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if ((max1 > 1.67597599124282407923e+153)){ + if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153))){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } + } - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; + if (double_tmp_result > eps){ + return false; } -#endif } else if(scz > bzmax) { double scz_bzmax = scz - bzmax; - if(max1 < scz_bzmax) + if(max1 < scz_bzmax){ max1 = scz_bzmax; + } distance += square(scz_bzmax); -#ifdef EARLY - { - double double_tmp_result = (distance - ssr); + double_tmp_result = (distance - ssr); - if ((max1 < 3.33558365626356687717e-147)){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if ((max1 > 1.67597599124282407923e+153)){ + if( ((max1 < 3.33558365626356687717e-147)) || ((max1 > 1.67597599124282407923e+153)) ){ + if(overestimate){ + return true; + }else{ CGAL_BRANCH_PROFILER_BRANCH_2(tmp); return Base::operator()(s,b); } - - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - - if (double_tmp_result > eps) - return false; } -#endif - } - double double_tmp_result = (distance - ssr); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); -#ifndef EARLY - if (max1 < 3.33558365626356687717e-147){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if (max1 > 1.67597599124282407923e+153){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); + if (double_tmp_result > eps){ + return false; + } } -#endif - double eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if (double_tmp_result > eps) - return false; - else - { - if (double_tmp_result < -eps) + // double_tmp_result and eps were growing all the time + // no need to test for > eps as done earlier in at least one case + if (double_tmp_result < -eps){ + return true; + } else { + if(overestimate){ return true; - else { - if(overestimate){ - return true; - } - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); } + CGAL_BRANCH_PROFILER_BRANCH_2(tmp); + return Base::operator()(s,b); } CGAL_BRANCH_PROFILER_BRANCH_2(tmp); From da17681a0630c01e2e1fa458117919c082bded89 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Wed, 10 Mar 2021 13:26:32 +0000 Subject: [PATCH 25/27] Avoid computing degree several times --- .../internal/Isotropic_remeshing/remesh_impl.h | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 551c1294cfaa..16f46e7cb7c4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -828,8 +828,6 @@ namespace internal { std::cout << "Equalize valences..." << std::endl; #endif -#define VALE -#ifdef VALE typedef typename boost::property_map >::type Vertex_degree; Vertex_degree degree = get(CGAL::dynamic_vertex_property_t(), mesh_); @@ -841,7 +839,6 @@ namespace internal { vertex_descriptor t = target(h, mesh_); put(degree, t, get(degree,t)+1); } -#endif unsigned int nb_flips = 0; for(edge_descriptor e : edges(mesh_)) @@ -855,17 +852,12 @@ namespace internal { vertex_descriptor vb = target(he, mesh_); vertex_descriptor vc = target(next(he, mesh_), mesh_); vertex_descriptor vd = target(next(opposite(he, mesh_), mesh_), mesh_); -#ifdef VALE + int vva = get(degree,va), tvva = target_valence(va); int vvb = get(degree, vb), tvvb = target_valence(vb); int vvc = get(degree,vc), tvvc = target_valence(vc); int vvd = get(degree,vd), tvvd = target_valence(vd); -#else - int vva = valence(va), tvva = target_valence(va); - int vvb = valence( vb), tvvb = target_valence(vb); - int vvc = valence(vc), tvvc = target_valence(vc); - int vvd = valence(vd), tvvd = target_valence(vd); -#endif + int deviation_pre = CGAL::abs(vva - tvva) + CGAL::abs(vvb - tvvb) + CGAL::abs(vvc - tvvc) @@ -1227,11 +1219,6 @@ namespace internal { out.close(); } - int valence(const vertex_descriptor& v) const - { - return static_cast(degree(v, mesh_)); - } - int target_valence(const vertex_descriptor& v) const { return (has_border_ && is_border(v, mesh_)) ? 4 : 6; From 109a893657c18d80216f6177f6301523f4557581 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 19 Mar 2021 11:30:05 +0100 Subject: [PATCH 26/27] remove #ifdef that should be defined --- .../internal/Isotropic_remeshing/remesh_impl.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 16f46e7cb7c4..fa93c69d8de5 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -875,12 +875,12 @@ namespace internal { vvb -= 1; vvc += 1; vvd += 1; -#ifdef VALE + put(degree, va, vva); put(degree, vb, vvb); put(degree, vc, vvc); put(degree, vd, vvd); -#endif + ++nb_flips; #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS @@ -916,17 +916,17 @@ namespace internal { CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) ); CGAL_assertion( !get(ecmap_, edge(he, mesh_)) ); CGAL::Euler::flip_edge(he, mesh_); -#ifdef VALE + vva += 1; vvb += 1; vvc -= 1; vvd -= 1; - put(degree, va, vva); - put(degree, vb, vvb); - put(degree, vc, vvc); - put(degree, vd, vvd); -#endif + put(degree, va, vva); + put(degree, vb, vvb); + put(degree, vc, vvc); + put(degree, vd, vvd); + --nb_flips; CGAL_assertion_code(Halfedge_status s3 = status(he)); From 29a40dc40b982465b42f72826590902d0f59aab8 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 19 Mar 2021 15:47:12 +0100 Subject: [PATCH 27/27] Remove the timer and the output of the result --- .../isotropic_remeshing_example.cpp | 24 ++----------------- 1 file changed, 2 insertions(+), 22 deletions(-) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp index bb246a9393c9..2fb29d976ef6 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_example.cpp @@ -1,7 +1,5 @@ -#define EARLY #include #include -#include #include #include #include @@ -37,9 +35,6 @@ int main(int argc, char* argv[]) { const char* filename = (argc > 1) ? argv[1] : "data/pig.off"; - CGAL::Timer t; - t.start(); - Mesh mesh; if(!PMP::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) { @@ -47,10 +42,7 @@ int main(int argc, char* argv[]) return 1; } - std::cout << "Read: " << t.time() << " sec" << std::endl; - t.reset(); - - double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 1.0; + double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 0.04; unsigned int nb_iter = 3; std::cout << "Split border..."; @@ -60,26 +52,14 @@ int main(int argc, char* argv[]) PMP::split_long_edges(border, target_edge_length, mesh); std::cout << "done." << std::endl; - - std::cout << "took " << t.time() << " sec" << std::endl; - t.reset(); - std::cout << "Start remeshing of " << filename << " (" << num_faces(mesh) << " faces)..." << std::endl; PMP::isotropic_remeshing(faces(mesh), target_edge_length, mesh, PMP::parameters::number_of_iterations(nb_iter) - .protect_constraints(true)); //i.e. protect border, here + .protect_constraints(true)); //i.e. protect border, here std::cout << "Remeshing done." << std::endl; - std::cout << "took " << t.time() << " sec" << std::endl; - t.reset(); - - std::ofstream out("out.off"); - out.precision(17); - //out << mesh << std::endl; - out.close(); - return 0; }