diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 613d08c174bf..7cabbdae1df7 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,6 +416,27 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } + 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() + (p, GeomTraits().compute_squared_distance_3_object()(p, bound)), bb,true)? + 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, Boolean_tag::value>()); + } + template CGAL::Comparison_result operator()(const Point& p, const Solid& pr, const FT& sq_distance) const { @@ -423,6 +446,7 @@ class AABB_traits CGAL::SMALLER : CGAL::LARGER; } + }; Closest_point closest_point_object() const {return Closest_point(*this);} diff --git a/BGL/include/CGAL/boost/graph/Euler_operations.h b/BGL/include/CGAL/boost/graph/Euler_operations.h index 8d5268da6e0d..d26d6a4bcebe 100644 --- a/BGL/include/CGAL/boost/graph/Euler_operations.h +++ b/BGL/include/CGAL/boost/graph/Euler_operations.h @@ -737,7 +737,7 @@ 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); 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..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 @@ -168,8 +168,9 @@ 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 { CGAL_BRANCH_PROFILER_3(std::string("semi-static failures/attempts/calls to : ") + std::string(CGAL_PRETTY_FUNCTION), tmp); @@ -177,30 +178,47 @@ 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); + 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; - + double double_tmp_result = 0; + double eps = 0; if(scx < bxmin) { double bxmin_scx = bxmin - scx; max1 = bxmin_scx; distance = square(bxmin_scx); + double_tmp_result = (distance - ssr); + + 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); + } + } + + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps){ + return false; + } } else if(scx > bxmax) { @@ -208,66 +226,136 @@ class Do_intersect_3 max1 = scx_bxmax; distance = square(scx_bxmax); + double_tmp_result = (distance - ssr); + + 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); + } + } + + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps){ + return false; + } } + if(scy < bymin) { double bymin_scy = bymin - scy; - if(max1 < bymin_scy) + if(max1 < bymin_scy){ max1 = bymin_scy; + } distance += square(bymin_scy); + double_tmp_result = (distance - ssr); + + 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); + } + } + + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps){ + return false; + } } else if(scy > bymax) { double scy_bymax = scy - bymax; - if(max1 < scy_bymax) + if(max1 < scy_bymax){ max1 = scy_bymax; - + } distance += square(scy_bymax); + double_tmp_result = (distance - ssr); + + 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); + } + } + + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps){ + return false; + } } + if(scz < bzmin) { double bzmin_scz = bzmin - scz; - if(max1 < bzmin_scz) + if(max1 < bzmin_scz){ max1 = bzmin_scz; - + } distance += square(bzmin_scz); + double_tmp_result = (distance - ssr); + + 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); + } + } + + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); + + if (double_tmp_result > eps){ + return false; + } } else if(scz > bzmax) { double scz_bzmax = scz - bzmax; - if(max1 < scz_bzmax) + if(max1 < scz_bzmax){ max1 = scz_bzmax; + } distance += square(scz_bzmax); - } + double_tmp_result = (distance - ssr); + + 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 double_tmp_result = (distance - ssr); - double max2 = CGAL::abs(ssr); + eps = 1.99986535548615598560e-15 * (std::max) (ssr, square(max1)); - if ((max1 < 3.33558365626356687717e-147) || (max2 < 1.11261183279326254436e-293)){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); - } - if ((max1 > 1.67597599124282407923e+153) || (max2 > 2.80889552322236673473e+306)){ - CGAL_BRANCH_PROFILER_BRANCH_2(tmp); - return Base::operator()(s,b); + if (double_tmp_result > eps){ + return false; + } } - double eps = 1.99986535548615598560e-15 * (std::max) (max2, (max1 * 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 { - 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); 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 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..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 @@ -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(certainly(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(certainly(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(certainly(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(certainly(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 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..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,6 +1,5 @@ #include #include - #include #include #include @@ -9,6 +8,7 @@ #include #include +#include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Surface_mesh Mesh; @@ -42,7 +42,7 @@ int main(int argc, char* argv[]) return 1; } - double target_edge_length = 0.04; + double target_edge_length = (argc > 2) ? std::stod(std::string(argv[2])) : 0.04; unsigned int nb_iter = 3; std::cout << "Split border..."; @@ -52,13 +52,12 @@ int main(int argc, char* argv[]) PMP::split_long_edges(border, target_edge_length, mesh); std::cout << "done." << std::endl; - 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; 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..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 @@ -30,6 +30,7 @@ #include #include #include +#include #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP # ifndef CGAL_PMP_COMPUTE_NORMAL_DEBUG @@ -90,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); @@ -106,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); } } @@ -506,6 +513,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()) @@ -669,7 +677,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 #include +#include #include #include #include @@ -43,6 +44,7 @@ #include #include #include +#include #include #include @@ -825,6 +827,19 @@ namespace internal { #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << "Equalize valences..." << std::endl; #endif + + 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); + } + unsigned int nb_flips = 0; for(edge_descriptor e : edges(mesh_)) { @@ -838,10 +853,11 @@ namespace internal { 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); + 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); + int deviation_pre = CGAL::abs(vva - tvva) + CGAL::abs(vvb - tvvb) + CGAL::abs(vvc - tvvc) @@ -859,6 +875,12 @@ namespace internal { vvb -= 1; vvc += 1; vvd += 1; + + put(degree, va, vva); + put(degree, vb, vvb); + put(degree, vc, vvc); + put(degree, vd, vvd); + ++nb_flips; #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS @@ -894,6 +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_); + + vva += 1; + vvb += 1; + vvc -= 1; + vvd -= 1; + + 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)); @@ -1186,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; @@ -1393,6 +1421,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 +1438,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; @@ -1879,8 +1913,10 @@ 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 initialized(nb_patches,false); + std::vector normal_per_patch(nb_patches); - std::vector< std::vector > normals_per_patch(nb_patches); for(halfedge_descriptor hd : hedges) { Halfedge_status s = status(hd); @@ -1892,18 +1928,18 @@ 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]){ + 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] = n; + initialized[index] = true; } return true; } @@ -1917,19 +1953,6 @@ namespace internal { 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_; 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) {