diff --git a/BGL/include/CGAL/boost/graph/Face_filtered_graph.h b/BGL/include/CGAL/boost/graph/Face_filtered_graph.h index 794329af4f2..25547b2ad4d 100644 --- a/BGL/include/CGAL/boost/graph/Face_filtered_graph.h +++ b/BGL/include/CGAL/boost/graph/Face_filtered_graph.h @@ -437,7 +437,7 @@ struct Face_filtered_graph initialize_halfedge_indices(); } - ///change the set of selected faces using a patch id + /// changes the set of selected faces using a patch id template void set_selected_faces(typename boost::property_traits::value_type face_patch_id, FacePatchIndexMap face_patch_index_map) @@ -466,7 +466,7 @@ struct Face_filtered_graph reset_indices(); } - /// change the set of selected faces using a range of patch ids + /// changes the set of selected faces using a range of patch ids template void set_selected_faces(const FacePatchIndexRange& selected_face_patch_indices, FacePatchIndexMap face_patch_index_map @@ -506,7 +506,7 @@ struct Face_filtered_graph reset_indices(); } - /// change the set of selected faces using a range of face descriptors + /// changes the set of selected faces using a range of face descriptors template void set_selected_faces(const FaceRange& selection) { @@ -569,18 +569,21 @@ struct Face_filtered_graph { return selected_halfedges[get(himap, halfedge(e,_graph))]; } - ///returns the number of selected faces - size_type number_of_faces()const + + /// returns the number of selected faces + size_type number_of_faces() const { return selected_faces.count(); } -///returns the number of selected vertices. - size_type number_of_vertices()const + + /// returns the number of selected vertices. + size_type number_of_vertices() const { return selected_vertices.count(); } -///returns the number of selected halfedges. - size_type number_of_halfedges()const + + /// returns the number of selected halfedges. + size_type number_of_halfedges() const { return selected_halfedges.count(); } @@ -612,18 +615,15 @@ struct Face_filtered_graph return bind_property_maps(himap, make_property_map(halfedge_indices) ); } - /// returns `true` if around any vertex of a selected face, - /// there is at most one connected set of selected faces. + /// returns `true` if around any vertex of a selected face there is at most a single umbrella bool is_selection_valid() const { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - // Non-manifoldness can appear either: - // - if 'pm' is pinched at a vertex. While traversing the incoming halfedges at this vertex, - // we will meet strictly more than one border halfedge. - // - if there are multiple umbrellas around a vertex. In that case, we will find a non-visited - // halfedge that has for target a vertex that is already visited. + // Non-manifoldness can appear if there are multiple umbrellas around a vertex. + // In that case, we will find a non-visited halfedge that has for target a vertex + // that is already visited. boost::unordered_set vertices_visited; boost::unordered_set halfedges_handled; @@ -643,15 +643,11 @@ struct Face_filtered_graph if(!vertices_visited.insert(vd).second) return false; - std::size_t border_halfedge_counter = 0; - - // Can't simply call halfedges_around_target(vd, *this) because 'halfedge(vd)' is not necessarily 'hd' + // Can't call halfedges_around_target(vd, *this) because 'halfedge(vd)' is not necessarily 'hd' halfedge_descriptor ihd = hd; do { halfedges_handled.insert(ihd); - if(is_border(ihd, *this)) - ++border_halfedge_counter; do { @@ -660,9 +656,6 @@ struct Face_filtered_graph while(!is_in_cc(ihd) && ihd != hd); } while(ihd != hd); - - if(border_halfedge_counter > 1) - return false; } return true; @@ -1050,7 +1043,8 @@ typename boost::graph_traits< Face_filtered_graph >: next(typename boost::graph_traits< Face_filtered_graph >::halfedge_descriptor h, const Face_filtered_graph & w) { - CGAL_assertion(w.is_in_cc(h)); + CGAL_precondition(w.is_in_cc(h)); + if(w.is_in_cc(next(h, w.graph()))) return next(h, w.graph()); @@ -1074,8 +1068,8 @@ typename boost::graph_traits< Face_filtered_graph >: prev(typename boost::graph_traits< Face_filtered_graph >::halfedge_descriptor h, const Face_filtered_graph & w) { + CGAL_precondition(w.is_in_cc(h)); - CGAL_assertion(w.is_in_cc(h)); if(w.is_in_cc(prev(h, w.graph()))) return prev(h, w.graph()); diff --git a/BGL/test/BGL/test_Face_filtered_graph.cpp b/BGL/test/BGL/test_Face_filtered_graph.cpp index 6b75dc7b1f9..099be0f4f69 100644 --- a/BGL/test/BGL/test_Face_filtered_graph.cpp +++ b/BGL/test/BGL/test_Face_filtered_graph.cpp @@ -484,8 +484,8 @@ void test_invalid_selections() face_range.push_back(SM::Face_index(2)); face_range.push_back(SM::Face_index(3)); - CGAL::Face_filtered_graph bad_fg(mesh, face_range); - assert(!bad_fg.is_selection_valid()); + CGAL::Face_filtered_graph pinched_fg(mesh, face_range); + assert(pinched_fg.is_selection_valid()); // this creates a non-manifold vertex (multiple umbrellas) clear(mesh); @@ -495,8 +495,8 @@ void test_invalid_selections() face_range.clear(); merge_vertices(SM::Vertex_index(1337), SM::Vertex_index(87), face_range, mesh); - CGAL::Face_filtered_graph bad_fg_2(mesh, face_range); - assert(!bad_fg_2.is_selection_valid()); + CGAL::Face_filtered_graph many_umbrellas_fg(mesh, face_range); + assert(!many_umbrellas_fg.is_selection_valid()); } int main() diff --git a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h index e95ccf8fe1c..f927d296338 100644 --- a/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h +++ b/Optimal_bounding_box/include/CGAL/Optimal_bounding_box/oriented_bounding_box.h @@ -92,7 +92,7 @@ void construct_oriented_bounding_box(const PointRange& points, obb_points[6] = cp(xmax, ymin, zmax); obb_points[7] = cp(xmax, ymax, zmax); - // Apply the inverse rotation to the rotated axis aligned bounding box + // Apply the inverse rotation to the rotated axis-aligned bounding box for(std::size_t i=0; i<8; ++i) { obb_points[i] = inverse_transformation.transform(obb_points[i]); @@ -141,13 +141,13 @@ void compute_best_transformation(const PointRange& points, rot(1, 0), rot(1, 1), rot(1, 2), rot(2, 0), rot(2, 1), rot(2, 2)); - // inverse transformation is simply the transposed since the matrix is unitary + // the inverse transformation is simply the transposed matrix since the matrix is unitary inverse_transformation = Aff_transformation_3(rot(0, 0), rot(1, 0), rot(2, 0), rot(0, 1), rot(1, 1), rot(2, 1), rot(0, 2), rot(1, 2), rot(2, 2)); } -// Following two functions are overloads to dispatch depending on return type +// The following two functions are overloads to dispatch depending on the return type template void construct_oriented_bounding_box(const PointRange& points, CGAL::Aff_transformation_3& transformation, @@ -325,6 +325,7 @@ void oriented_bounding_box(const PointRange& points, { using CGAL::parameters::choose_parameter; using CGAL::parameters::get_parameter; + using CGAL::parameters::is_default_parameter; typedef typename CGAL::GetPointMap::type PointMap; @@ -350,7 +351,8 @@ void oriented_bounding_box(const PointRange& points, const unsigned int seed = choose_parameter(get_parameter(np, internal_np::random_seed), -1); // undocumented CGAL::Random fixed_seed_rng(seed); - CGAL::Random& rng = (seed == unsigned(-1)) ? CGAL::get_default_random() : fixed_seed_rng; + CGAL::Random& rng = is_default_parameter(get_parameter(np, internal_np::random_seed)) ? + CGAL::get_default_random() : fixed_seed_rng; #ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG std::cout << "Random seed: " << rng.get_seed() << std::endl; diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h index f4eb6980549..bf369321af4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polygon_mesh.h @@ -46,16 +46,21 @@ struct Tracer_polyhedron int i, int k, bool last = true) { - if(i + 1 == k) { return P[i+1]; } + if(i + 1 == k) + return P[i+1]; halfedge_descriptor h, g; - if(i+2 == k){ + if(i+2 == k) + { if(last) - { - h = P[i+1]; - Euler::fill_hole(h,pmesh); } + { + h = P[i + 1]; + Euler::fill_hole(h, pmesh); + } else - { h = Euler::add_face_to_border(prev(P[i+1],pmesh), P[i+2/*k*/], pmesh); } + { + h = Euler::add_face_to_border(prev(P[i + 1], pmesh), P[i + 2 /*k*/], pmesh); + } CGAL_assertion(face(h,pmesh) != boost::graph_traits::null_face()); *out++ = face(h,pmesh); @@ -68,12 +73,14 @@ struct Tracer_polyhedron g = operator()(lambda, la, k, false); if(last) - { - h = g; - Euler::fill_hole(g,pmesh); - } + { + h = g; + Euler::fill_hole(g, pmesh); + } else - { h = Euler::add_face_to_border(prev(h,pmesh), g, pmesh); } + { + h = Euler::add_face_to_border(prev(h, pmesh), g, pmesh); + } CGAL_assertion(face(h,pmesh) != boost::graph_traits::null_face()); *out++ = face(h,pmesh); @@ -106,28 +113,29 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh, typedef std::map Vertex_map; typedef typename Vertex_map::iterator Vertex_map_it; - #ifdef CGAL_PMP_HOLE_FILLING_DEBUG +#ifdef CGAL_PMP_HOLE_FILLING_DEBUG CGAL::Timer timer; timer.start(); - #endif +#endif - std::vector P, Q; + std::vector P, Q; std::vector P_edges; Vertex_map vertex_map; int id = 0; Hedge_around_face_circulator circ(border_halfedge,pmesh), done(circ); - do{ + do + { P.push_back(get(vpmap, target(*circ, pmesh))); Q.push_back(get(vpmap, target(next(opposite(next(*circ,pmesh),pmesh),pmesh),pmesh))); P_edges.push_back(*circ); - if(!vertex_map.insert(std::make_pair(target(*circ,pmesh), id++)).second) { - #ifndef CGAL_TEST_SUITE + if(!vertex_map.insert(std::make_pair(target(*circ,pmesh), id++)).second) + { +#ifndef CGAL_TEST_SUITE CGAL_warning_msg(false, "Returning no output. Non-manifold vertex is found on boundary!"); - #else +#else std::cerr << "W: Returning no output. Non-manifold vertex is found on boundary!\n"; - #endif - return std::make_pair(out, - CGAL::internal::Weight_min_max_dihedral_and_area::NOT_VALID()); +#endif + return std::make_pair(out, CGAL::internal::Weight_min_max_dihedral_and_area::NOT_VALID()); } } while (++circ != done); @@ -151,52 +159,52 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh, if(v_it_neigh_it != vertex_map.end()) //other endpoint found in the map { int v_it_neigh_id = v_it_neigh_it->second; - if( v_it_neigh_id != v_it_prev && v_it_neigh_id != v_it_next ) - { //there is an edge incident to v_it, which is not next or previous + if(v_it_neigh_id != v_it_prev && v_it_neigh_id != v_it_next) + { + //there is an edge incident to v_it, which is not next or previous //from vertex_map (checked by comparing IDs) if(v_it_id < v_it_neigh_id) // to include each edge only once - { existing_edges.push_back(std::make_pair(v_it_id, v_it_neigh_id)); } + existing_edges.push_back(std::make_pair(v_it_id, v_it_neigh_id)); } } } while(++circ_vertex != done_vertex); } - //#define CGAL_USE_WEIGHT_INCOMPLETE - #ifdef CGAL_USE_WEIGHT_INCOMPLETE +//#define CGAL_USE_WEIGHT_INCOMPLETE +#ifdef CGAL_USE_WEIGHT_INCOMPLETE typedef CGAL::internal::Weight_calculator, CGAL::internal::Is_valid_existing_edges_and_degenerate_triangle> WC; - #else +#else typedef CGAL::internal::Weight_calculator WC; - #endif +#endif CGAL::internal::Is_valid_existing_edges_and_degenerate_triangle is_valid(existing_edges); // fill hole using polyline function, with custom tracer for PolygonMesh - Tracer_polyhedron - tracer(out, pmesh, P_edges); + Tracer_polyhedron tracer(out, pmesh, P_edges); #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2 - if(use_cdt && triangulate_hole_polyline_with_cdt(P, tracer, is_valid, k, max_squared_distance)) - { - return std::make_pair(tracer.out, CGAL::internal::Weight_min_max_dihedral_and_area(0,0)); - } + if(use_cdt && triangulate_hole_polyline_with_cdt(P, tracer, is_valid, k, max_squared_distance)) + return std::make_pair(tracer.out, CGAL::internal::Weight_min_max_dihedral_and_area(0,0)); #endif CGAL::internal::Weight_min_max_dihedral_and_area weight = - triangulate_hole_polyline(P, Q, tracer, WC(is_valid), - use_delaunay_triangulation, k) -#ifdef CGAL_USE_WEIGHT_INCOMPLETE - .weight // get actual weight in Weight_incomplete +#ifndef CGAL_USE_WEIGHT_INCOMPLETE + triangulate_hole_polyline(P, Q, tracer, WC(is_valid), use_delaunay_triangulation, k); +#else + // get actual weight in Weight_incomplete + triangulate_hole_polyline(P, Q, tracer, WC(is_valid), use_delaunay_triangulation, k).weight; #endif - ; - #ifdef CGAL_PMP_HOLE_FILLING_DEBUG +#ifdef CGAL_PMP_HOLE_FILLING_DEBUG std::cerr << "Hole filling: " << timer.time() << " sc." << std::endl; timer.reset(); - #endif +#endif + return std::make_pair(tracer.out, weight); } -}// namespace internal -}// namespace Polygon_mesh_processing -}// namespace CGAL +} // namespace internal +} // namespace Polygon_mesh_processing +} // namespace CGAL + #endif //CGAL_HOLE_FILLING_TRIANGULATE_HOLE_POLYHEDRON_3_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h index 5eab37058cd..403494666ce 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Triangulate_hole_polyline.h @@ -1492,11 +1492,11 @@ triangulate_hole_polyline(const PointRange1& points, typedef Kernel K; typedef typename K::Point_3 Point_3; - #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 +#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 typedef CGAL::internal::Triangulate_hole_polyline_DT Fill_DT; - #else +#else CGAL_USE(use_delaunay_triangulation); - #endif +#endif typedef CGAL::internal::Triangulate_hole_polyline Fill; std::vector P(boost::begin(points), boost::end(points)); @@ -1510,20 +1510,19 @@ triangulate_hole_polyline(const PointRange1& points, } typename WeightCalculator::Weight w = - #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 +#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 use_delaunay_triangulation ? Fill_DT().operator()(P,Q,tracer,WC) : - #endif +#endif Fill().operator()(P,Q,tracer,WC); #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 - if (use_delaunay_triangulation - && w == WeightCalculator::Weight::NOT_VALID()) + if(use_delaunay_triangulation && w == WeightCalculator::Weight::NOT_VALID()) w = Fill().operator()(P, Q, tracer, WC); #endif - #ifdef CGAL_PMP_HOLE_FILLING_DEBUG +#ifdef CGAL_PMP_HOLE_FILLING_DEBUG std::cerr << w << std::endl; - #endif +#endif return w; } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine.h index 03bdf356a37..effb3fb7d3f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine.h @@ -27,7 +27,7 @@ namespace Polygon_mesh_processing { /*! \ingroup PMP_meshing_grp - @brief refines a region of a triangle mesh + @brief refines a region of a triangle mesh. @tparam TriangleMesh model of `MutableFaceGraph` @tparam FaceRange range of face descriptors, model of `Range`. diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h index 6d34f67bd9a..030864e6bbe 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h @@ -105,7 +105,7 @@ namespace Polygon_mesh_processing { * * \cgalParamNBegin{protect_constraints} * \cgalParamDescription{If `true`, the edges set as constrained in `edge_is_constrained_map` -* (or by default the boundary edges) are not split nor collapsed during remeshing.} +* (or by default the boundary edges) are neither split nor collapsed during remeshing.} * \cgalParamType{Boolean} * \cgalParamDefault{`false`} * \cgalParamExtra{Note that around constrained edges that have their length higher than diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h index 18b4b07940e..938b1b64f5a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h @@ -15,10 +15,12 @@ #include +#include #include #include #include #include +#include #include #include #include @@ -26,6 +28,9 @@ #include #endif +#include +#include +#include #include #include #include @@ -33,14 +38,17 @@ #include #include #include +#ifdef CGAL_PMP_REPAIR_SI_USE_OBB_IN_COMPACTIFICATION +#include +#endif #include - #include #include #include #include #include #include +#include #include #include #include @@ -75,44 +83,45 @@ static int self_intersections_solved_by_unconstrained_hole_filling = 0; // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- -// @todo these could be extracted to somewhere else, it's useful in itself -template -FaceOutputIterator replace_faces_with_patch(const std::vector::vertex_descriptor>& border_vertices, - const std::set::vertex_descriptor>& interior_vertices, - const std::vector::halfedge_descriptor>& border_hedges, - const std::set::edge_descriptor>& interior_edges, - const std::set::face_descriptor>& faces, - const std::vector >& patch, - PolygonMesh& pmesh, - VPM& vpm, - FaceOutputIterator out) +template +FaceOutputIterator replace_faces_with_patch_without_reuse(const std::vector::vertex_descriptor>& border_vertices, + const std::set::face_descriptor>& faces, + const std::vector >& patch, + PolygonMesh& pmesh, + VertexPointMap vpm, + FaceOutputIterator out) { - CGAL_static_assertion((std::is_same::value_type, Point>::value)); - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::edge_descriptor edge_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef std::vector Point_face; typedef std::vector Vertex_face; - CGAL_precondition(is_valid_polygon_mesh(pmesh)); + std::map point_to_vd; - // To be used to create new elements - std::vector vertex_stack(interior_vertices.begin(), interior_vertices.end()); - std::vector edge_stack(interior_edges.begin(), interior_edges.end()); - std::vector face_stack(faces.begin(), faces.end()); - - // Introduce new vertices, convert the patch in vertex patches - std::vector patch_with_vertices; - patch_with_vertices.reserve(patch.size()); - - std::map point_to_vs; - - // first, add those for which the vertex will not change + // First, add those for which the vertex will not change for(const vertex_descriptor v : border_vertices) - point_to_vs[get(vpm, v)] = v; + { + // In this version, remove_face() will get rid of isolated vertices so only vertices incident + // to at least one face that is not going to be removed will be kept + bool kept_vertex = false; + for(face_descriptor f : faces_around_target(halfedge(v, pmesh), pmesh)) + { + if(f != boost::graph_traits::null_face() && faces.count(f) == 0) + { + kept_vertex = true; + break; + } + } + + if(kept_vertex) + point_to_vd[get(vpm, v)] = v; + } + + for(face_descriptor f : faces) + Euler::remove_face(halfedge(f, pmesh), pmesh); + + CGAL_assertion(is_valid_polygon_mesh(pmesh)); // now build a correspondence map and the faces with vertices const vertex_descriptor null_v = boost::graph_traits::null_vertex(); @@ -125,10 +134,95 @@ FaceOutputIterator replace_faces_with_patch(const std::vector::iterator it; - std::tie(it, success) = point_to_vs.insert(std::make_pair(p, null_v)); + std::tie(it, success) = point_to_vd.emplace(p, null_v); vertex_descriptor& v = it->second; - if(success) // first time we meet that point, means it`s an interior point and we need to make a new vertex + if(success) + { + // first time we meet that point, means it's an interior point and we need to make a new vertex + v = add_vertex(pmesh); + put(vpm, v, p); + } + + vface.push_back(v); + } + + face_descriptor new_f = boost::graph_traits::null_face(); + if(Euler::can_add_face(vface, pmesh)) + new_f = Euler::add_face(vface, pmesh); + + if(new_f == boost::graph_traits::null_face()) + { + std::cerr << "Error: failed to insert patch face??" << std::endl; + return out; + } + + out++ = new_f; + } + + return out; +} + + +// @todo these could be extracted to somewhere else, it's useful in itself +template +FaceOutputIterator replace_faces_with_patch(const std::vector::vertex_descriptor>& border_vertices, + const std::set::vertex_descriptor>& interior_vertices, + const std::vector::halfedge_descriptor>& border_hedges, + const std::set::edge_descriptor>& interior_edges, + const std::set::face_descriptor>& faces, + const std::vector >& patch, + PolygonMesh& pmesh, + VertexPointMap vpm, + FaceOutputIterator out) +{ + CGAL_static_assertion((std::is_same::value_type, Point>::value)); + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef std::vector Point_face; + typedef std::vector Vertex_face; + + CGAL_precondition(is_valid_polygon_mesh(pmesh)); + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Replacing range with patch: "; + std::cout << faces.size() << " triangles removed, " << patch.size() << " created\n"; +#endif + + // To be used to create new elements + std::vector vertex_stack(interior_vertices.begin(), interior_vertices.end()); + std::vector edge_stack(interior_edges.begin(), interior_edges.end()); + std::vector face_stack(faces.begin(), faces.end()); + + // Introduce new vertices, convert the patch in vertex patches + std::vector patch_with_vertices; + patch_with_vertices.reserve(patch.size()); + + std::map point_to_vd; + + // first, add those for which the vertex will not change + for(const vertex_descriptor v : border_vertices) + point_to_vd[get(vpm, v)] = v; + + // now build a correspondence map and the faces with vertices + const vertex_descriptor null_v = boost::graph_traits::null_vertex(); + for(const Point_face& face : patch) + { + Vertex_face vface; + vface.reserve(face.size()); + + for(const Point& p : face) + { + bool success; + typename std::map::iterator it; + std::tie(it, success) = point_to_vd.emplace(p, null_v); + vertex_descriptor& v = it->second; + + if(success) // first time we meet that point, means it's an interior point and we need to make a new vertex { if(vertex_stack.empty()) { @@ -160,16 +254,13 @@ FaceOutputIterator replace_faces_with_patch(const std::vector::null_face(); -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::vector new_faces; -#endif for(const Vertex_face& vface : patch_with_vertices) { @@ -183,10 +274,8 @@ FaceOutputIterator replace_faces_with_patch(const std::vector::null_face()); *out++ = f; -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - new_faces.push_back(f); -#endif std::vector hedges; hedges.reserve(vface.size()); @@ -199,8 +288,8 @@ FaceOutputIterator replace_faces_with_patch(const std::vector::null_halfedge())); + std::tie(it, success) = halfedge_map.emplace(std::make_pair(vi, vj), + boost::graph_traits::null_halfedge()); halfedge_descriptor& h = it->second; if(success) // this halfedge is an interior halfedge @@ -244,29 +333,20 @@ FaceOutputIterator replace_faces_with_patch(const std::vector +template FaceOutputIterator replace_faces_with_patch(const std::set::face_descriptor>& face_range, const std::vector >& patch, PolygonMesh& pmesh, - VPM& vpm, + VertexPointMap vpm, FaceOutputIterator out) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -317,11 +397,11 @@ FaceOutputIterator replace_faces_with_patch(const std::set +template void replace_faces_with_patch(const std::set::face_descriptor>& faces, const std::vector >& patch, PolygonMesh& pmesh, - VPM& vpm) + VertexPointMap vpm) { CGAL::Emptyset_iterator out; replace_faces_with_patch(faces, patch, pmesh, vpm, out); @@ -329,39 +409,17 @@ void replace_faces_with_patch(const std::set -void back_up_face_range_as_point_patch(std::vector >& point_patch, - const FaceRange& face_range, - const PolygonMesh& tmesh, - const VertexPointMap vpm) -{ - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - point_patch.reserve(face_range.size()); - - for(const face_descriptor f : face_range) - { - std::vector face_points; - for(const halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh)) - face_points.push_back(get(vpm, target(h, tmesh))); - - point_patch.push_back(face_points); - } -} - -// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- - template -void constrain_sharp_and_border_edges(const FaceRange& faces, - TriangleMesh& tmesh, - EdgeConstrainMap& eif, - const bool constrain_sharp_edges, - const double dihedral_angle, - const double /*weak_DA*/, - VertexPointMap vpm, - const GeomTraits& gt) +void constrain_edges(const FaceRange& faces, + TriangleMesh& tmesh, + const bool constrain_border_edges, + const bool constrain_sharp_edges, + const double dihedral_angle, + const double /*weak_DA*/, + EdgeConstrainMap& eif, + VertexPointMap vpm, + const GeomTraits& gt) { typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; @@ -370,16 +428,18 @@ void constrain_sharp_and_border_edges(const FaceRange& faces, typedef typename GeomTraits::FT FT; typedef typename GeomTraits::Vector_3 Vector; - std::map is_border_of_selection; + std::unordered_map is_border_of_selection; for(face_descriptor f : faces) { - // @fixme what about nm vertices for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh)) { // Default initialization is guaranteed to be `false`. Thus, meet it once will switch // the value to `true` and meeting it twice will switch back to `false`. const edge_descriptor e = edge(h, tmesh); - is_border_of_selection[e] = !(is_border_of_selection[e]); + if(constrain_sharp_edges) + is_border_of_selection[e] = !(is_border_of_selection[e]); + else + is_border_of_selection[e] = false; } } @@ -387,10 +447,7 @@ void constrain_sharp_and_border_edges(const FaceRange& faces, CGAL::Polygon_mesh_processing::experimental::detect_sharp_edges_pp(faces, tmesh, dihedral_angle, eif, parameters::weak_dihedral_angle(weak_DA)); - // borders are also constrained - for(const auto& ep : is_border_of_selection) - if(ep.second) - put(eif, ep.first, true); + // ... #else // this is basically the code that is in detect_features (at the very bottom) // but we do not want a folding to be marked as a sharp feature so the dihedral angle is also @@ -401,6 +458,9 @@ void constrain_sharp_and_border_edges(const FaceRange& faces, for(const auto& ep : is_border_of_selection) { bool flag = ep.second; + if(!constrain_border_edges) + flag = false; + if(constrain_sharp_edges && !flag) { const halfedge_descriptor h = halfedge(ep.first, tmesh); @@ -409,7 +469,7 @@ void constrain_sharp_and_border_edges(const FaceRange& faces, const face_descriptor f1 = face(h, tmesh); const face_descriptor f2 = face(opposite(h, tmesh), tmesh); - // @todo cache normals + // @speed cache normals const Vector n1 = compute_face_normal(f1, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)); const Vector n2 = compute_face_normal(f2, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)); const FT c = gt.compute_scalar_product_3_object()(n1, n2); @@ -477,8 +537,8 @@ bool remove_self_intersections_with_smoothing(std::set > patch; for(const face_descriptor f : faces(local_mesh)) { @@ -515,7 +576,6 @@ bool remove_self_intersections_with_smoothing(std::set new_faces; replace_faces_with_patch(face_range, patch, tmesh, vpm, std::inserter(new_faces, new_faces.end())); - CGAL_assertion(!does_self_intersect(new_faces, tmesh, parameters::vertex_point_map(vpm))); #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG if(constrain_sharp_edges) @@ -561,61 +621,33 @@ bool order_border_halfedge_range(std::vector -void dump_cc(const FaceContainer& cc_faces, - const TriangleMesh& mesh, - const std::string filename) -{ - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename GetVertexPointMap::const_type VertexPointMap; - VertexPointMap vpm =get_const_property_map(vertex_point, mesh); - - std::ofstream out(filename); - out.precision(17); - - out << "OFF\n"; - out << 3*cc_faces.size() << " " << cc_faces.size() << " 0\n"; - - for(const face_descriptor f : cc_faces) - { - out << get(vpm, source(halfedge(f, mesh), mesh)) << "\n"; - out << get(vpm, target(halfedge(f, mesh), mesh)) << "\n"; - out << get(vpm, target(next(halfedge(f, mesh), mesh), mesh)) << "\n"; - } - - int id = 0; - for(const face_descriptor f : cc_faces) - { - CGAL_USE(f); - out << "3 " << id << " " << id+1 << " " << id+2 << "\n"; - id += 3; - } - - out.close(); -} - template -void dump_tentative_hole(std::vector >& point_patch, - const std::string filename) +void dump_patch(const std::string filename, + std::vector >& point_patch) { std::ofstream out(filename); out << std::setprecision(17); +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Writing " << point_patch.size() << " face(s) into " << filename << std::endl; +#endif + + std::vector points; + std::vector > faces; + std::map unique_points_with_id; for(const std::vector& face : point_patch) for(const Point& p : face) - unique_points_with_id.insert(std::make_pair(p, 0)); + unique_points_with_id.emplace(p, 0); out << "OFF\n"; out << unique_points_with_id.size() << " " << point_patch.size() << " 0\n"; int unique_id = 0; - for(auto& pp : unique_points_with_id) + for(auto& e : unique_points_with_id) { - out << pp.first << "\n"; - pp.second = unique_id++; + e.second = unique_id++; + out << e.first << "\n"; } for(const std::vector& face : point_patch) @@ -630,6 +662,48 @@ void dump_tentative_hole(std::vector >& point_patch, out.close(); } +template +void dump_cc(const std::string filename, + const FaceContainer& cc_faces, + const PolygonMesh& pmesh, + const VertexPointMap vpm) +{ + typedef typename boost::property_traits::value_type Point; + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Writing " << cc_faces.size() << " face(s) into " << filename << std::endl; +#endif + + std::unordered_map vertex_ids; + std::stringstream vss, fss; + std::size_t id = 0; + for(face_descriptor f : cc_faces) + { + fss << degree(f, pmesh); + for(vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh)) + { + auto res = vertex_ids.emplace(v, id); + if(res.second) // insert was successful (first time seeing this vertex) + { + ++id; + vss << get(vpm, v) << "\n"; + } + + fss << " " << res.first->second /*id*/; + } + fss << "\n"; + } + + std::ofstream out(filename); + out << std::setprecision(17); + + out << "OFF\n"; + out << id << " " << cc_faces.size() << " 0\n"; + out << vss.str() << "\n" << fss.str() << std::endl; +} #endif // CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- @@ -677,118 +751,55 @@ construct_artificial_third_point(const typename boost::graph_traits -bool construct_tentative_hole_patch(std::vector::vertex_descriptor>& cc_border_vertices, - std::set::vertex_descriptor>& cc_interior_vertices, - std::set::edge_descriptor>& cc_interior_edges, - const std::vector& hole_points, - const std::vector& third_points, - const std::vector::halfedge_descriptor>& cc_border_hedges, - const std::set::face_descriptor>& cc_faces, - std::vector >& point_patch, - const TriangleMesh& tmesh, - VertexPointMap /*vpm*/, - const GeomTraits& /*gt*/) +template +bool check_patch_compatibility(const std::vector >& patch, + const std::vector::vertex_descriptor>& border_vertices, + const std::vector::halfedge_descriptor>& border_hedges, + const std::set::edge_descriptor>& interior_edges, + const TriangleMesh& tmesh, + const VertexPointMap vpm) { - CGAL_static_assertion((std::is_same::value_type, Point>::value)); + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; + std::map point_to_vd; + for(vertex_descriptor v : border_vertices) + point_to_vd[get(vpm, v)] = v; - typedef CGAL::Triple Face_indices; - - CGAL_assertion(cc_border_hedges.size() == cc_border_hedges.size()); - CGAL_assertion(hole_points.size() == third_points.size()); - - // Collect vertices and edges inside the current selection cc: first collect all vertices and - // edges incident to the faces to remove... - for(const face_descriptor f : cc_faces) - { - for(halfedge_descriptor h : halfedges_around_face(halfedge(f, tmesh), tmesh)) - { - if(halfedge(target(h, tmesh), tmesh) == h) // limit the number of insertions - cc_interior_vertices.insert(target(h, tmesh)); - - cc_interior_edges.insert(edge(h, tmesh)); - } - } - - // ... and then remove those on the boundary - for(halfedge_descriptor h : cc_border_hedges) - { - cc_interior_vertices.erase(target(h, tmesh)); - cc_interior_edges.erase(edge(h, tmesh)); - } - - // try to triangulate the hole using default parameters - // (using Delaunay search space if CGAL_HOLE_FILLING_DO_NOT_USE_DT3 is not defined) - std::vector hole_faces; - if(hole_points.size() > 3) - triangulate_hole_polyline(hole_points, third_points, std::back_inserter(hole_faces)); - else - hole_faces.emplace_back(0, 1, 2); // trivial hole filling - - if(hole_faces.empty()) - { -#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 - #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: Failed to fill a hole using Delaunay search space.\n"; - #endif - - triangulate_hole_polyline(hole_points, third_points, std::back_inserter(hole_faces), - parameters::use_delaunay_triangulation(false)); -#endif - if(hole_faces.empty()) - { -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: Failed to fill a hole using the whole search space.\n"; -#endif - return false; - } - } - -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT - std::cout << " DEBUG: " << hole_faces.size() << " faces in the patch" << std::endl; - std::vector > to_dump; - for(const Face_indices& face : hole_faces) - { - to_dump.emplace_back(std::initializer_list{hole_points[face.first], - hole_points[face.second], - hole_points[face.third]}); - } - - CGAL_assertion(to_dump.size() == hole_faces.size()); - - static int hole_id = 0; - std::stringstream oss; - oss << "results/tentative_hole_" << hole_id++ << ".off" << std::ends; - const std::string filename = oss.str().c_str(); - - dump_tentative_hole(to_dump, filename); -#endif - - // make sure that the hole filling is valid, we check that no - // edge already in the mesh is present in hole_faces. + // make sure that the hole filling is valid: check that no edge + // already in the mesh is present in hole_faces. bool non_manifold_edge_found = false; - for(const Face_indices& triangle : hole_faces) + for(const std::vector& f : patch) { - std::array edges = make_array(triangle.first, triangle.second, - triangle.second, triangle.third, - triangle.third, triangle.first); - for(int k=0; k<3; ++k) + for(int i=0; i<3; ++i) { - const int vi = edges[2*k], vj = edges[2*k+1]; + const Point& p0 = f[i]; + const Point& p1 = f[(i+1)%3]; - // ignore boundary edges - if(vi+1 == vj || (vj == 0 && static_cast(vi) == cc_border_vertices.size()-1)) + auto p0_it = point_to_vd.find(p0); + auto p1_it = point_to_vd.find(p1); + + // @fixme + // If any of the vertices is an inner point created through refine(), we don't have an easy way + //to know the possible correspondency with an existing vertex of the mesh: it might be a vertex + // part of a completely different CC. Unfortunately, a nm edge could be created with this vertex, + // but the complexity to check all vertices of the mesh is horrible (even spatially filtered, + // this needs to be updated, ...) + if(p0_it == point_to_vd.end() || p1_it == point_to_vd.end()) continue; - halfedge_descriptor h = halfedge(cc_border_vertices[vi], cc_border_vertices[vj], tmesh).first; - if(h != boost::graph_traits::null_halfedge() && - cc_interior_edges.count(edge(h, tmesh)) == 0) + const vertex_descriptor v0 = p0_it->second; + const vertex_descriptor v1 = p1_it->second; + + halfedge_descriptor h = halfedge(v0, v1, tmesh).first; // null halfedge if not found + if(h != boost::graph_traits::null_halfedge()) { - non_manifold_edge_found = true; - break; + if(std::find(border_hedges.begin(), border_hedges.end(), h) == border_hedges.end() && + interior_edges.count(edge(h, tmesh)) == 0) + { + non_manifold_edge_found = true; + break; + } } } @@ -805,33 +816,379 @@ bool construct_tentative_hole_patch(std::vector +bool check_patch_sanity(const std::vector >& patch) +{ + std::set > unique_faces; + std::map, int> unique_edges; + + for(const std::vector& face : patch) { - point_patch.emplace_back(std::initializer_list{hole_points[face.first], - hole_points[face.second], - hole_points[face.third]}); + if(!unique_faces.emplace(face.begin(), face.end()).second) // this face had already been found + return false; + + int i = (unique_edges.insert(std::make_pair(std::set { face[0], face[1] }, 0)).first->second)++; + if(i == 2) // non-manifold edge + return false; + + i = (unique_edges.insert(std::make_pair(std::set { face[1], face[2] }, 0)).first->second)++; + if(i == 2) // non-manifold edge + return false; + + i = (unique_edges.insert(std::make_pair(std::set { face[2], face[0] }, 0)).first->second)++; + if(i == 2) // non-manifold edge + return false; } + // Check for self-intersections within the patch + // @todo something better than just making a mesh out of the soup? + std::vector points; + std::vector > faces; + std::map ids; + + std::size_t c = 0; + for(const std::vector& face : patch) + { + std::vector ps_f; + for(const Point& pt : face) + { + std::size_t id = c; + auto is_insert_successful = ids.emplace(pt, c); + if(is_insert_successful.second) // first time we've seen that point + { + ++c; + points.push_back(pt); + } + else // already seen that point + { + id = is_insert_successful.first->second; + } + + CGAL_assertion(id < points.size()); + ps_f.push_back(id); + } + + faces.push_back(ps_f); + } + + TriangleMesh patch_mesh; + if(is_polygon_soup_a_polygon_mesh(faces)) + polygon_soup_to_polygon_mesh(points, faces, patch_mesh); + else + return false; + + if(does_self_intersect(patch_mesh)) + { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: Found acceptable hole-filling patch.\n"; + std::cout << " DEBUG: Tentative patch has self-intersections." << std::endl; +#endif + + return false; + } + + return true; +} + +template +bool construct_hole_patch(std::vector >& hole_faces, + const std::vector& hole_points, + const std::vector& third_points, + const GeomTraits& gt) +{ + if(hole_points.size() > 3) + { + triangulate_hole_polyline(hole_points, third_points, std::back_inserter(hole_faces), + parameters::geom_traits(gt)); + } + else + { + hole_faces.emplace_back(0, 1, 2); // trivial hole filling + } + + if(hole_faces.empty()) + { +#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 + #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to fill a hole using Delaunay search space.\n"; + #endif + + triangulate_hole_polyline(hole_points, third_points, std::back_inserter(hole_faces), + parameters::use_delaunay_triangulation(false).geom_traits(gt)); +#endif + if(hole_faces.empty()) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to fill a hole using the whole search space.\n"; +#endif + return false; + } + } + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT + std::cout << " DEBUG: " << hole_faces.size() << " faces in the patch" << std::endl; + std::vector > to_dump; + for(const auto& face : hole_faces) + { + to_dump.emplace_back(std::initializer_list{hole_points[face.first], + hole_points[face.second], + hole_points[face.third]}); + } + + CGAL_assertion(to_dump.size() == hole_faces.size()); + + static int patch_id = 0; + std::stringstream oss; + oss << "results/raw_patch_" << patch_id++ << ".off" << std::ends; + const std::string filename = oss.str().c_str(); + + dump_patch(filename, to_dump); #endif return true; } +template +struct Mesh_projection_functor +{ + typedef typename GeomTraits::Point_3 Point_3; + typedef typename GeomTraits::Triangle_3 Triangle_3; + + typedef std::vector Triangle_container; + typedef CGAL::AABB_triangle_primitive Primitive; + typedef CGAL::AABB_traits Traits; + typedef CGAL::AABB_tree Tree; + + template + Mesh_projection_functor(const TriangleMesh& mesh, + const VPM vpm) + { + triangles.reserve(num_faces(mesh)); + for(auto f : faces(mesh)) + triangles.emplace_back(get(vpm, target(halfedge(f, mesh), mesh)), + get(vpm, target(next(halfedge(f, mesh), mesh), mesh)), + get(vpm, source(halfedge(f, mesh), mesh))); + + tree.insert(std::cbegin(triangles), std::cend(triangles)); + } + + Point_3 operator()(const Point_3& p) const { return tree.closest_point(p); } + +private: + Triangle_container triangles; + Tree tree; +}; + +template +bool adapt_patch(std::vector >& point_patch, + const Projector& projector, + const TriangleMesh&, + const GeomTraits&) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GeomTraits::FT FT; + + std::vector soup_points; + std::vector > soup_faces; + + std::size_t pid = 0; + std::map point_ids; + for(const auto& fp : point_patch) + { + CGAL_assertion(fp.size() == 3); + std::array f; + for(std::size_t i=0; i<3; ++i) + { + auto res = point_ids.emplace(fp[i], pid); + if(res.second) + { + soup_points.push_back(fp[i]); + ++pid; + } + f[i] = res.first->second; + } + soup_faces.push_back(f); + } + + TriangleMesh local_mesh; + auto local_vpm = get(vertex_point, local_mesh); + + polygon_soup_to_polygon_mesh(soup_points, soup_faces, local_mesh); + bool has_SI = does_self_intersect(local_mesh); + + std::vector border_hedges; + border_halfedges(faces(local_mesh), local_mesh, std::back_inserter(border_hedges)); + typename TriangleMesh::template Property_map selected_edge = // @fixme BGL + local_mesh.template add_property_map("e:selected",false).first; + + for(halfedge_descriptor h : border_hedges) + selected_edge[edge(h, local_mesh)] = true; + + std::vector new_vertices; + refine(local_mesh, faces(local_mesh), CGAL::Emptyset_iterator(), std::back_inserter(new_vertices)); + + for(vertex_descriptor v : new_vertices) + local_mesh.point(v) = projector(get(local_vpm, v)); // @fixme BGL + + // The projector can create degenerate faces + if(!remove_degenerate_faces(local_mesh)) + return !has_SI; + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT + static int adapted_patch_id = 0; + std::stringstream oss; + oss << "results/adapted_patch_" << adapted_patch_id++ << ".off" << std::ends; + const std::string filename = oss.str().c_str(); + std::cout << " DEBUG: Writing " << point_patch.size() << " faces into " << filename << std::endl; + IO::write_polygon_mesh(filename, local_mesh); +#endif + + // If the adapted tentative patch has SI, revert back to the base patch + if(does_self_intersect(local_mesh)) + return !has_SI; // if the base patch also has self-intersections, we are done + + // Replace the tentative patch with the new, self-intersection-less, adapted patch + point_patch.clear(); + point_patch.reserve(num_faces(local_mesh)); + + for(face_descriptor f : faces(local_mesh)) + { + std::vector fp { get(local_vpm, target(halfedge(f, local_mesh), local_mesh)), + get(local_vpm, target(next(halfedge(f, local_mesh), local_mesh), local_mesh)), + get(local_vpm, source(halfedge(f, local_mesh), local_mesh)) }; + point_patch.push_back(fp); + } + + return true; +} + +// This overload uses hole filling to construct a patch and tests the manifoldness of the patch +template +bool construct_manifold_hole_patch(std::vector >& point_patch, + const std::vector& hole_points, + const std::vector& third_points, + const std::vector::vertex_descriptor>& cc_border_vertices, + const std::vector::halfedge_descriptor>& cc_border_hedges, + const std::set::edge_descriptor>& cc_interior_edges, + const Projector& projector, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + typedef CGAL::Triple Face_indices; + + // Try to triangulate the hole using default parameters + // (using Delaunay search space if CGAL_HOLE_FILLING_DO_NOT_USE_DT3 is not defined) + std::vector hole_faces; + construct_hole_patch(hole_faces, hole_points, third_points, gt); + + std::vector > local_point_patch; + local_point_patch.reserve(hole_faces.size()); + for(const Face_indices& face : hole_faces) + { + local_point_patch.emplace_back(std::initializer_list{hole_points[face.first], + hole_points[face.second], + hole_points[face.third]}); + } + + if(!adapt_patch(local_point_patch, projector, tmesh, gt)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to adapt the patch..." << std::endl; +#endif + return false; + } + + // Check manifoldness compatibility with the rest of the mesh + if(!check_patch_compatibility(local_point_patch, cc_border_vertices, cc_border_hedges, cc_interior_edges, tmesh, vpm)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Incompatible patch" << std::endl; +#endif + return false; + } + + point_patch.reserve(point_patch.size() + local_point_patch.size()); + std::move(std::begin(local_point_patch), std::end(local_point_patch), std::back_inserter(point_patch)); + + bool is_sane = check_patch_sanity(point_patch); +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + if(is_sane) + std::cout << " DEBUG: Found sane hole-filling patch (" << point_patch.size() << " faces)\n"; + else + std::cout << " DEBUG: Insane hole-filling patch\n"; +#endif + + return is_sane; +} + +// This overloads fill the containers `cc_interior_vertices` and `cc_interior_edges` +template +bool construct_tentative_hole_patch_with_border(std::vector >& point_patch, + const std::vector& hole_points, + const std::vector& third_points, + const std::vector::vertex_descriptor>& cc_border_vertices, + const std::vector::halfedge_descriptor>& cc_border_hedges, + std::set::vertex_descriptor>& cc_interior_vertices, + std::set::edge_descriptor>& cc_interior_edges, + const std::set::face_descriptor>& cc_faces, + const Projector& projector, + const TriangleMesh& tmesh, + const VertexPointMap vpm, + const GeomTraits& gt) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + CGAL_assertion(hole_points.size() == third_points.size()); + + // Collect vertices and edges inside the current selection cc: first collect all vertices and + // edges incident to the faces to remove... + for(const face_descriptor f : cc_faces) + { + for(halfedge_descriptor h : halfedges_around_face(halfedge(f, tmesh), tmesh)) + { + if(halfedge(target(h, tmesh), tmesh) == h) // to limit the number of insertions + cc_interior_vertices.insert(target(h, tmesh)); + + cc_interior_edges.insert(edge(h, tmesh)); + } + } + + // ... and then remove those on the boundary + for(halfedge_descriptor h : cc_border_hedges) + { + cc_interior_vertices.erase(target(h, tmesh)); + cc_interior_edges.erase(edge(h, tmesh)); + } + + return construct_manifold_hole_patch(point_patch, hole_points, third_points, + cc_border_vertices, cc_border_hedges, cc_interior_edges, + projector, tmesh, vpm, gt); +} + // This function constructs the ranges `hole_points` and `third_points`. Note that for a sub-hole, // these two ranges are constructed in another function because we don't want to set 'third_points' // for edges that are on the border of the sub-hole but not on the border of the (full) hole. -template -bool construct_tentative_hole_patch(std::vector::vertex_descriptor>& cc_border_vertices, +template +bool construct_tentative_hole_patch(std::vector::value_type> >& patch, + std::vector::vertex_descriptor>& cc_border_vertices, std::set::vertex_descriptor>& cc_interior_vertices, std::set::edge_descriptor>& cc_interior_edges, const std::vector::halfedge_descriptor>& cc_border_hedges, const std::set::face_descriptor>& cc_faces, - std::vector::value_type> >& patch, + const Projector& projector, const TriangleMesh& tmesh, - VertexPointMap vpm, + const VertexPointMap vpm, const GeomTraits& gt) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -839,6 +1196,8 @@ bool construct_tentative_hole_patch(std::vector::value_type Point; + cc_border_vertices.reserve(cc_border_hedges.size()); + std::vector hole_points, third_points; hole_points.reserve(cc_border_hedges.size()); third_points.reserve(cc_border_hedges.size()); @@ -859,19 +1218,21 @@ bool construct_tentative_hole_patch(std::vector= 3); - return construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges, - hole_points, third_points, cc_border_hedges, cc_faces, - patch, tmesh, vpm, gt); + return construct_tentative_hole_patch_with_border(patch, hole_points, third_points, + cc_border_vertices, cc_border_hedges, + cc_interior_vertices, cc_interior_edges, + cc_faces, projector, tmesh, vpm, gt); } -// In that overload, we don't know the border of the patch because the face range is a sub-region +// In this overload, we don't know the border of the patch because the face range is a sub-region // of the hole. We also construct `hole_points` and `third_points`, but with no third point for internal // sharp edges because a local self-intersection is usually caused by folding and thus we do not want -// a third point resulting from folding to constrain the way we fill the hole in the wrong way. -template +// a third point resulting from folding to wrongly influence the hole filling process. +template bool construct_tentative_sub_hole_patch(std::vector::value_type> >& patch, const std::set::face_descriptor>& sub_cc_faces, const std::set::face_descriptor>& cc_faces, + const Projector& projector, TriangleMesh& tmesh, VertexPointMap vpm, const GeomTraits& gt) @@ -921,8 +1282,8 @@ bool construct_tentative_sub_hole_patch(std::vector cc_interior_vertices; std::set cc_interior_edges; @@ -953,99 +1314,25 @@ bool construct_tentative_sub_hole_patch(std::vector -bool check_patch_sanity(const std::vector >& patch) -{ - std::set > unique_faces; - std::map, int> unique_edges; - - for(const std::vector& face : patch) - { - if(!unique_faces.emplace(face.begin(), face.end()).second) // this face had already been found - return false; - - int i = (unique_edges.insert(std::make_pair(std::set { face[0], face[1] }, 0)).first->second)++; - if(i == 2) // non-manifold edge - return false; - - i = (unique_edges.insert(std::make_pair(std::set { face[1], face[2] }, 0)).first->second)++; - if(i == 2) // non-manifold edge - return false; - - i = (unique_edges.insert(std::make_pair(std::set { face[2], face[0] }, 0)).first->second)++; - if(i == 2) // non-manifold edge - return false; - } - - // Check for self-intersections - // Don't know anything better than just making a mesh out of the soup for now... - std::vector points; - std::vector > faces; - std::map ids; - - std::size_t c = 0; - for(const std::vector& face : patch) - { - std::vector ps_f; - for(const Point& pt : face) - { - std::size_t id = c; - std::pair::iterator, bool> is_insert_successful = - ids.insert(std::make_pair(pt, c)); - if(is_insert_successful.second) // first time we've seen that point - { - ++c; - points.push_back(pt); - } - else // already seen that point - { - id = is_insert_successful.first->second; - } - - CGAL_assertion(id < points.size()); - ps_f.push_back(id); - } - - faces.push_back(ps_f); - } - - TriangleMesh patch_mesh; - if(is_polygon_soup_a_polygon_mesh(faces)) - polygon_soup_to_polygon_mesh(points, faces, patch_mesh); - else - return false; - - if(does_self_intersect(patch_mesh)) - { -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: Tentative patch has self-intersections." << std::endl; -#endif - - return false; - } - - return true; -} - // This function is only called when the hole is NOT subdivided into smaller holes -template +template bool fill_hole(std::vector::halfedge_descriptor>& cc_border_hedges, - std::set::face_descriptor>& cc_faces, + const std::set::face_descriptor>& cc_faces, std::set::face_descriptor>& working_face_range, - TriangleMesh& tmesh, const PolyhedralEnvelope& cc_envelope, + const Projector& projector, + TriangleMesh& tmesh, VertexPointMap vpm, - const GeomTraits& gt) + const GeomTraits& gt, + bool reuse_faces = true) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; @@ -1059,6 +1346,10 @@ bool fill_hole(std::vector::halfedge_ if(!order_border_halfedge_range(cc_border_hedges, tmesh)) { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to orient the boundary??\n"; +#endif + CGAL_assertion(false); // we shouldn't fail to orient the boundary cycle of the complete hole return false; } @@ -1070,9 +1361,8 @@ bool fill_hole(std::vector::halfedge_ cc_border_vertices.reserve(cc_border_hedges.size()); std::vector > patch; - if(!construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges, - cc_border_hedges, cc_faces, patch, tmesh, vpm, gt) || - !check_patch_sanity(patch)) + if(!construct_tentative_hole_patch(patch, cc_border_vertices, cc_interior_vertices, cc_interior_edges, + cc_border_hedges, cc_faces, projector, tmesh, vpm, gt)) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << " DEBUG: Failed to find acceptable hole patch\n"; @@ -1081,7 +1371,7 @@ bool fill_hole(std::vector::halfedge_ return false; } - if (!cc_envelope.is_empty() && !cc_envelope(patch)) + if(!cc_envelope.is_empty() && !cc_envelope(patch)) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << " DEBUG: Patch is not inside the input polyhedral envelope\n"; @@ -1089,22 +1379,32 @@ bool fill_hole(std::vector::halfedge_ return false; } - // Could renew the range directly within the patch replacement function - // to avoid erasing and re-adding the same face +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Found acceptable hole-filling patch (" << patch.size() << " faces)\n"; +#endif + for(const face_descriptor f : cc_faces) working_face_range.erase(f); // Plug the new triangles in the mesh, reusing previous edges and faces - replace_faces_with_patch(cc_border_vertices, cc_interior_vertices, - cc_border_hedges, cc_interior_edges, - cc_faces, patch, tmesh, vpm, - std::inserter(working_face_range, working_face_range.end())); + if(reuse_faces) + { + replace_faces_with_patch(cc_border_vertices, cc_interior_vertices, + cc_border_hedges, cc_interior_edges, + cc_faces, patch, tmesh, vpm, + std::inserter(working_face_range, working_face_range.end())); + } + else + { + replace_faces_with_patch_without_reuse(cc_border_vertices, cc_faces, patch, tmesh, vpm, + std::inserter(working_face_range, working_face_range.end())); + } -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT_INTERMEDIATE_FULL_MESH static int filed_hole_id = 0; std::stringstream oss; oss << "results/filled_basic_" << filed_hole_id++ << ".off" << std::ends; - std::ofstream(oss.str().c_str()) << std::setprecision(17) << tmesh; + CGAL::IO::write_polygon_mesh(oss.str().c_str(), tmesh, CGAL::parameters::stream_precision(17)); #endif CGAL_postcondition(is_valid_polygon_mesh(tmesh)); @@ -1113,11 +1413,12 @@ bool fill_hole(std::vector::halfedge_ } // Same function as above but border of the hole is not known -template -bool fill_hole(std::set::face_descriptor>& cc_faces, +template +bool fill_hole(const std::set::face_descriptor>& cc_faces, std::set::face_descriptor>& working_face_range, - TriangleMesh& tmesh, const PolyhedralEnvelope& cc_envelope, + const Projector& projector, + TriangleMesh& tmesh, VertexPointMap vpm, const GeomTraits& gt) { @@ -1138,20 +1439,20 @@ bool fill_hole(std::set::face_descrip } if(order_border_halfedge_range(cc_border_hedges, tmesh)) - return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, - cc_envelope,vpm, gt); + return fill_hole(cc_border_hedges, cc_faces, working_face_range, cc_envelope, projector, tmesh, vpm, gt); else return false; } -template +template bool fill_hole_with_constraints(std::vector::halfedge_descriptor>& cc_border_hedges, - std::set::face_descriptor>& cc_faces, + const std::set::face_descriptor>& cc_faces, std::set::face_descriptor>& working_face_range, TriangleMesh& tmesh, const double dihedral_angle, const double weak_DA, const PolyhedralEnvelope& cc_envelope, + const Projector& projector, VertexPointMap vpm, const GeomTraits& gt) { @@ -1168,7 +1469,8 @@ bool fill_hole_with_constraints(std::vector::type EIFMap; EIFMap eif = get(Edge_property_tag(), tmesh); - constrain_sharp_and_border_edges(cc_faces, tmesh, eif, true /*constrain_sharp_edges*/, dihedral_angle, weak_DA, vpm, gt); + constrain_edges(cc_faces, tmesh, true /*constrain_border_edges*/, true /*constrain_sharp_edges*/, + dihedral_angle, weak_DA, eif, vpm, gt); // Partition the hole using these constrained edges std::set visited_faces; @@ -1180,86 +1482,54 @@ bool fill_hole_with_constraints(std::vector sub_cc; Polygon_mesh_processing::connected_component(f, tmesh, std::inserter(sub_cc, sub_cc.end()), CGAL::parameters::edge_is_constrained_map(eif)); visited_faces.insert(sub_cc.begin(), sub_cc.end()); -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << "CC of size " << sub_cc.size() << " (total: " << cc_faces.size() << ")" << std::endl; -#endif ++cc_counter; #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT - dump_cc(sub_cc, tmesh, "results/current_cc.off"); + dump_cc("results/current_cc.off", sub_cc, tmesh, vpm); #endif // The mesh is not modified, but 'patch' gets filled - if(!construct_tentative_sub_hole_patch(patch, sub_cc, cc_faces, tmesh, vpm, gt)) + if(!construct_tentative_sub_hole_patch(patch, sub_cc, cc_faces, projector, tmesh, vpm, gt)) { // Something went wrong while finding a potential cover for the a sub-hole --> use basic hole-filling - return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, - cc_envelope,vpm, gt); + return fill_hole(cc_border_hedges, cc_faces, working_face_range, cc_envelope, projector, tmesh, vpm, gt); } } -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << cc_counter << " independent sub holes" << std::endl; -#endif + #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT - std::ofstream out("results/hole_fillers.off"); - out.precision(17); - - out << "OFF\n"; - out << 3*patch.size() << " " << patch.size() << " 0\n"; - - for(const auto& f : patch) - { - for(const auto& pt : f) - out << pt << "\n"; - } - - int id = 0; - for(std::size_t i=0; i(patch)) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << "Unhealthy patch, use base fill_hole" << std::endl; + std::cout << " DEBUG: Unhealthy patch, defaulting to basic fill_hole" << std::endl; #endif - return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, - cc_envelope, vpm, gt); + return fill_hole(cc_border_hedges, cc_faces, working_face_range, cc_envelope, projector, tmesh, vpm, gt); } // check if the patch is inside the input polyhedral envelope if(!cc_envelope.is_empty() && !cc_envelope(patch)) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << "Patch is not entirely inside the input polyhedral envelope, use base fill_hole" << std::endl; + std::cout << " DEBUG: Patch is not entirely inside the input polyhedral envelope, defaulting to basic fill_hole" << std::endl; #endif - return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, - cc_envelope, vpm, gt); + return fill_hole(cc_border_hedges, cc_faces, working_face_range, cc_envelope, projector, tmesh, vpm, gt); } - // Plug the hole-filling patch in the mesh - std::set new_faces; - replace_faces_with_patch(cc_faces, patch, tmesh, vpm, std::inserter(new_faces, new_faces.end())); - - // Otherwise it should have failed the sanity check - CGAL_assertion(!does_self_intersect(new_faces, tmesh, parameters::vertex_point_map(vpm))); - - // Update working range with the new faces for(const face_descriptor f : cc_faces) working_face_range.erase(f); - working_face_range.insert(new_faces.begin(), new_faces.end()); + // Plug the hole-filling patch in the mesh + replace_faces_with_patch(cc_faces, patch, tmesh, vpm, + std::inserter(working_face_range, working_face_range.end())); return true; } @@ -1352,15 +1622,15 @@ bool is_simple_3(const std::vector::h return true; } -template +template bool remove_self_intersections_with_hole_filling(std::vector::halfedge_descriptor>& cc_border_hedges, - std::set::face_descriptor>& cc_faces, + const std::set::face_descriptor>& cc_faces, std::set::face_descriptor>& working_face_range, TriangleMesh& tmesh, - bool local_self_intersection_removal, const double strong_dihedral_angle, const double weak_dihedral_angle, const PolyhedralEnvelope& cc_envelope, + const Projector& projector, VertexPointMap vpm, const GeomTraits& gt) { @@ -1380,69 +1650,284 @@ bool remove_self_intersections_with_hole_filling(std::vector +template +bool handle_CC_with_complex_topology(std::vector::halfedge_descriptor>& cc_border_hedges, + const std::set::face_descriptor>& cc_faces, + std::set::face_descriptor>& working_face_range, + TriangleMesh& tmesh, + const double strong_dihedral_angle, + const double weak_dihedral_angle, + const bool preserve_genus, + const PolyhedralEnvelope& cc_envelope, + const Projector& projector, + VertexPointMap vpm, + const GeomTraits& gt) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GeomTraits::FT FT; + typedef typename boost::property_traits::value_type Point; + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: CC with Euler_chi != 1" << std::endl; +#endif + + if(preserve_genus) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: CC not handled, selection is not a topological disk (preserve_genus=true)\n"; +#endif + return false; + } + + const CGAL::Face_filtered_graph ccmesh(tmesh, cc_faces); + if(!ccmesh.is_selection_valid()) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: invalid FFG selection\n"; +#endif + return false; + } + + std::vector boundary_reps; + extract_boundary_cycles(ccmesh, std::back_inserter(boundary_reps)); + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: " << boundary_reps.size() << " borders in the CC\n"; +#endif + + if(boundary_reps.size() == 1) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Complex topology but single border --> standard hole filling\n"; +#endif + + // If there is a single border, fill the hole as if it were a topological disk. + // This will lose some information since chi != -1, but preserve_genus = false here + return remove_self_intersections_with_hole_filling(cc_border_hedges, cc_faces, working_face_range, + tmesh, strong_dihedral_angle, weak_dihedral_angle, + cc_envelope, projector, vpm, gt); + } + + // From there on, there is more than one border + + std::vector is_hole_incident_to_patch(boundary_reps.size()); + std::vector hole_lengths(boundary_reps.size()); + + int holes_incident_to_patches_n = 0; + for(std::size_t hole_id = 0; hole_id longest_border_hedges; + halfedge_descriptor bh = boundary_reps[longest_border_id], end = bh; + do + { + longest_border_hedges.push_back(opposite(bh, tmesh)); + bh = prev(bh, ccmesh); // prev because we insert the opposite + } + while(bh != end); + + // 'false' because we can't do on-the-fly patching due to multiple boundary cycles + // @todo this currently doesn't attempt to constrain sharp edges + return fill_hole(longest_border_hedges, cc_faces, working_face_range, cc_envelope, projector, + tmesh, vpm, gt, false /*reuse*/); + } + + // If there exists some boundary cycles with "fake" border halfedges, hole-fill those +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Complex topology, some fake borders (" << holes_incident_to_patches_n << ")\n"; +#endif + + // This is needed for the patch insertion process at the end + std::vector all_border_vertices; // border vertices for all the borders to be filled + + // The patch is built iteratively and made of as many CCs as there are holes being filled + std::vector > patch; + + for(std::size_t hole_id=0; hole_id border_hedges; + halfedge_descriptor bh = boundary_reps[hole_id], end = bh; + do + { + border_hedges.push_back(opposite(bh, tmesh)); + bh = prev(bh, ccmesh); // prev because we insert the opposite + } + while(bh != end); + + std::vector border_vertices; + border_vertices.reserve(border_hedges.size()); + all_border_vertices.reserve(all_border_vertices.size() + border_hedges.size()); + + std::vector hole_points, third_points; + hole_points.reserve(border_hedges.size()); + third_points.reserve(border_hedges.size()); + + for(const halfedge_descriptor h : border_hedges) + { + CGAL_assertion(!is_border(h, tmesh)); + + const vertex_descriptor v = source(h, tmesh); + hole_points.push_back(get(vpm, v)); + + border_vertices.push_back(v); + all_border_vertices.push_back(v); + + if(is_border_edge(h, tmesh)) // h is incident to a real face + third_points.push_back(construct_artificial_third_point(h, tmesh, vpm, gt)); + else + third_points.push_back(get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh))); + } + + std::set interior_vertices; + std::set interior_edges; + + if(!construct_tentative_hole_patch_with_border(patch, hole_points, third_points, + border_vertices, border_hedges, + interior_vertices, interior_edges, + cc_faces, projector, tmesh, vpm, gt)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to fill hole #" << hole_id << "\n"; +#endif + + return false; + } + } + + // Built the patch from all the boundary cycles, put it in +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT + dump_patch("results/multiple_real_borders.off", patch); +#endif + + // We're assembling multiple patches so we could have the same face appearing multiple times... + if(!check_patch_sanity(patch)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Unhealthy patch(s), defaulting to basic fill_hole" << std::endl; +#endif + return false; + } + + // check if the patch is inside the input polyhedral envelope + if(!cc_envelope.is_empty() && !cc_envelope(patch)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Patch is not entirely inside the input polyhedral envelope, defaulting to basic fill_hole" << std::endl; +#endif + return false; + } + + for(const face_descriptor f : cc_faces) + working_face_range.erase(f); + + // Plug the hole-filling patch in the mesh + replace_faces_with_patch_without_reuse(all_border_vertices, cc_faces, patch, tmesh, vpm, + std::inserter(working_face_range, working_face_range.end())); + + return true; +} + +// the parameter `step` controls how many extra layers of faces we take around the range `faces_to_treat` +template std::pair -remove_self_intersections_one_step(std::set::face_descriptor>& faces_to_remove, +remove_self_intersections_one_step(std::set::face_descriptor>& faces_to_treat, std::set::face_descriptor>& working_face_range, TriangleMesh& tmesh, const int step, const bool preserve_genus, - const bool only_treat_self_intersections_locally, + const bool treat_all_CCs, const double strong_dihedral_angle, const double weak_dihedral_angle, const double containment_epsilon, + const Projector& projector, VertexPointMap vpm, const GeomTraits& gt, Visitor& visitor) { - typedef boost::graph_traits graph_traits; - typedef typename graph_traits::vertex_descriptor vertex_descriptor; - typedef typename graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename graph_traits::face_descriptor face_descriptor; + typedef boost::graph_traits graph_traits; + typedef typename graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename graph_traits::face_descriptor face_descriptor; - std::set faces_to_remove_copy = faces_to_remove; - -#if defined(CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG) || defined(CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT) - static int call_id = -1; - ++call_id; +#ifdef CGAL_PMP_REPAIR_SI_USE_OBB_IN_COMPACTIFICATION + typedef typename graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::property_traits::value_type Point; #endif -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << "##### running remove_self_intersections_one_step (#" << call_id << "), step " << step - << " with " << faces_to_remove.size() << " intersecting faces\n"; -#endif + std::set faces_to_treat_copy = faces_to_treat; bool something_was_done = false; // indicates if a region was successfully remeshed bool all_fixed = true; // indicates if all removal went well @@ -1451,34 +1936,41 @@ remove_self_intersections_one_step(std::set cc_faces; - std::vector queue(1, *faces_to_remove.begin()); // temporary queue + std::vector queue(1, *faces_to_treat.begin()); // temporary queue cc_faces.insert(queue.back()); while(!queue.empty()) { @@ -1488,9 +1980,9 @@ remove_self_intersections_one_step(std::set::null_face()) + if(adjacent_face != boost::graph_traits::null_face()) { - if(faces_to_remove.count(adjacent_face) != 0 && cc_faces.insert(adjacent_face).second) + if(faces_to_treat.count(adjacent_face) != 0 && cc_faces.insert(adjacent_face).second) queue.push_back(adjacent_face); } @@ -1499,21 +1991,17 @@ remove_self_intersections_one_step(std::set stack_for_expension; + +#ifdef CGAL_PMP_REPAIR_SI_USE_OBB_IN_COMPACTIFICATION + std::set cc_points; + for(face_descriptor f : cc_faces) + for(vertex_descriptor v : vertices_around_face(halfedge(f, tmesh), tmesh)) + cc_points.insert(get(vpm, v)); + + typedef typename GeomTraits::Aff_transformation_3 Aff_transformation; + Aff_transformation tr{CGAL::Identity_transformation()}; + + if(cc_points.size() > 3) + CGAL::oriented_bounding_box(cc_points, tr, CGAL::parameters::random_seed(0)); + + // Construct the rotated OBB Bbox_3 bb; + for(const Point& p : cc_points) + bb += (tr.transform(p)).bbox(); + +#else + Bbox_3 bb; +#endif + for(face_descriptor fd : cc_faces) { for(halfedge_descriptor h : halfedges_around_face(halfedge(fd, tmesh), tmesh)) { +#ifndef CGAL_PMP_REPAIR_SI_USE_OBB_IN_COMPACTIFICATION bb += get(vpm, target(h, tmesh)).bbox(); +#endif face_descriptor nf = face(opposite(h, tmesh), tmesh); if(nf != boost::graph_traits::null_face() && cc_faces.count(nf) == 0) - { stack_for_expension.push_back(opposite(h, tmesh)); - } } } @@ -1554,7 +2065,11 @@ remove_self_intersections_one_step(std::set > is_selected(cc_faces); + expand_face_selection_for_removal(cc_faces, tmesh, is_selected); #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << " DEBUG: " << cc_faces.size() << " faces in expanded and compactified CC\n"; #endif - // remove faces from the set to process - for(const face_descriptor f : cc_faces) - faces_to_remove.erase(f); - #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT - fname="results/expanded_compactified_r_"+std::to_string(call_id)+"_cc_"+std::to_string(cc_id)+"_s_"+std::to_string(step)+".off"; - std::cout << " DEBUG: Writing expanded and compactified CC #" << cc_id << " in " << fname << std::endl; - dump_cc(cc_faces, tmesh, fname); + fname = "results/expanded_compactified_step_"+std::to_string(step)+"_CC_"+std::to_string(cc_id)+".off"; + dump_cc(fname, cc_faces, tmesh, vpm); #endif - //Check for non-manifold vertices in the selection and remove them by selecting all incident faces: - // extract the set of halfedges that is on the boundary of the holes to be - // made. In addition, we make sure no hole to be created contains a vertex - // visited more than once along a hole border (pinched surface) - // We save the size of boundary_hedges to make sur halfedges added - // from non_filled_hole are not removed. - bool non_manifold_vertex_remaining_in_selection = false; - do - { - bool non_manifold_vertex_removed = false; //here non-manifold is for the 1D polyline - std::vector boundary_hedges; - for(face_descriptor fh : cc_faces) - { - halfedge_descriptor h = halfedge(fh, tmesh); - for(int i=0; i<3; ++i) - { - if(is_border(opposite(h, tmesh), tmesh) || cc_faces.count(face(opposite(h, tmesh), tmesh)) == 0) - boundary_hedges.push_back(h); + // Now, we have a proper selection to work on. - h = next(h, tmesh); - } - } + for(const face_descriptor f : cc_faces) + faces_to_treat.erase(f); - // detect vertices visited more than once along - // a hole border. We then remove all faces incident - // to such a vertex to force the removal of the vertex. - // Actually even if two holes are sharing a vertex, this - // vertex will be removed. It is not needed but since - // we do not yet have one halfedge per hole it is simpler - // and does not harm - std::set border_vertices; - for(halfedge_descriptor h : boundary_hedges) - { - if(!border_vertices.insert(target(h, tmesh)).second) - { - bool any_face_added = false; - for(halfedge_descriptor hh : halfedges_around_target(h, tmesh)) - { - if(!is_border(hh, tmesh)) - { - // add the face to the current selection - any_face_added |= cc_faces.insert(face(hh, tmesh)).second; - faces_to_remove.erase(face(hh, tmesh)); - } - } - - if(any_face_added) - non_manifold_vertex_removed = true; - else - non_manifold_vertex_remaining_in_selection = true; - } - } - - if(!non_manifold_vertex_removed) - break; - } - while(true); - - if(preserve_genus && non_manifold_vertex_remaining_in_selection) - { - topology_issue = true; -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: CC not handled due to the presence at least one non-manifold vertex\n"; - - ++unsolved_self_intersections; -#endif - - visitor.end_component_handling(); - continue; // cannot replace a patch containing a nm vertex by a disk - } - - // before running this function if preserve_genus=false, we duplicated - // all of them - CGAL_assertion(!non_manifold_vertex_remaining_in_selection); - - // Collect halfedges on the boundary of the region to be selected - // (pointing inside the domain to be remeshed) - std::vector cc_border_hedges; - for(face_descriptor fd : cc_faces) - { - halfedge_descriptor h = halfedge(fd, tmesh); - for(int i=0; i<3; ++i) - { - if(is_border(opposite(h, tmesh), tmesh) || cc_faces.count(face(opposite(h, tmesh), tmesh)) == 0) - cc_border_hedges.push_back(h); - - h = next(h, tmesh); - } - } - - if(cc_faces.size() == 1) // it is a triangle nothing better can be done + if(cc_faces.size() == 1) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - ++unsolved_self_intersections; + std::cout << " DEBUG: Compactified CC of size 1, moving on\n"; #endif visitor.end_component_handling(); continue; } - working_face_range.insert(cc_faces.begin(), cc_faces.end()); + bool self_intersects = does_self_intersect(cc_faces, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)); +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + if(!self_intersects) + std::cout << " DEBUG: No self-intersection within the CC\n"; +#endif - // Now, we have a proper selection that we can work on. + if(!treat_all_CCs && !self_intersects) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + ++unsolved_self_intersections; +#endif + + all_fixed = false; + visitor.end_component_handling(); + continue; + } #ifndef CGAL_PMP_REMOVE_SELF_INTERSECTION_NO_POLYHEDRAL_ENVELOPE_CHECK Polyhedral_envelope cc_envelope; - if (containment_epsilon!=0) + if(containment_epsilon != 0) cc_envelope = Polyhedral_envelope(cc_faces, tmesh, containment_epsilon); #else - struct Return_true { + struct Return_true + { bool is_empty() const { return true; } bool operator()(const std::vector >&) const { return true; } bool operator()(const TriangleMesh&) const { return true; } @@ -1718,188 +2145,120 @@ remove_self_intersections_one_step(std::set mesh_non_border_hedges; - std::set mesh_border_hedge; - - for(halfedge_descriptor h : cc_border_hedges) - { - if(!is_border(opposite(h, tmesh), tmesh)) - mesh_non_border_hedges.push_back(h); - else - mesh_border_hedge.insert(opposite(h, tmesh)); - } - - if (mesh_border_hedge.empty()) - { -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: CC not handled, selection is not a topological disk (preserve_genus=false)\n"; - ++unsolved_self_intersections; -#endif - topology_issue = true; - continue; - } - - // we look for cycles of border halfedges and update selection_chi - while(!mesh_border_hedge.empty()) - { - // we must count the number of cycle of boundary edges - halfedge_descriptor h_b = *mesh_border_hedge.begin(), h=h_b; - mesh_border_hedge.erase(mesh_border_hedge.begin()); - do - { - h = next(h, tmesh); - if(h == h_b) - { - // found a cycle - selection_chi += 1; - break; - } - else - { - typename std::set::iterator it = mesh_border_hedge.find(h); - if(it == mesh_border_hedge.end()) - break; // not a cycle does not count - - mesh_border_hedge.erase(it); - } - } - while(true); - } - - if(selection_chi!=1) - { -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: CC not handled, selection is not a topological disk even if" - << " boundary cycles are removed: chi=" << selection_chi << "\n"; - ++unsolved_self_intersections; -#endif - topology_issue = true; + something_was_done = true; visitor.end_component_handling(); continue; } + #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG else { - cc_border_hedges.swap(mesh_non_border_hedges); + std::cout << " DEBUG: Could not be solved via smoothing\n"; + } + #endif + } - // count the number of cycles of halfedges of the boundary - std::map bhs; - for(halfedge_descriptor h : cc_border_hedges) - bhs[source(h, tmesh)] = target(h, tmesh); +#endif // ndef CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_SMOOTHING - int nbc=0; - while(!bhs.empty()) - { - ++nbc; - std::pair top=*bhs.begin(); - bhs.erase(bhs.begin()); - - do - { - typename std::map::iterator it_find = bhs.find(top.second); - if(it_find == bhs.end()) - break; - - top = *it_find; - bhs.erase(it_find); - } - while(true); - } - - if(nbc != 1) - { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: CC not handled because it is not a topological disk(" - << nbc << " boundary cycles)\n"; - ++unsolved_self_intersections; + std::cout << " DEBUG: Trying hole-filling based approach...\n"; #endif - all_fixed = false; - visitor.end_component_handling(); - continue; - } - else - { -#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << " DEBUG: CC that is not a topological disk but has only one boundary cycle(preserve_genus=false)\n"; -#endif - } + // Collect halfedges on the boundary of the region to be selected + // (incident to faces that are part of the CC) + std::vector cc_border_hedges; + for(face_descriptor fd : cc_faces) + { + for(halfedge_descriptor h : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + if(is_border(opposite(h, tmesh), tmesh) || cc_faces.count(face(opposite(h, tmesh), tmesh)) == 0) + cc_border_hedges.push_back(h); } } + // Whichever step we are at, no border means no expansion will change this selection + // This CC was not fixed by smoothing, and there is nothing hole filling can do + // @todo just remove the CC? + if(cc_border_hedges.empty()) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: CC is closed!\n"; // @todo wrap? + ++unsolved_self_intersections; +#endif + + all_fixed = false; + visitor.end_component_handling(); + continue; + } + + int selection_chi = euler_characteristic_of_selection(cc_faces, tmesh); + if(selection_chi != 1) // not a topological disk + { + if(!handle_CC_with_complex_topology(cc_border_hedges, cc_faces, working_face_range, + tmesh, strong_dihedral_angle, weak_dihedral_angle, + preserve_genus, cc_envelope, projector, vpm, gt)) + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: Failed to handle complex CC\n"; + ++unsolved_self_intersections; +#endif + topology_issue = true; + all_fixed = false; + } + else + { +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + ++self_intersections_solved_by_unconstrained_hole_filling; +#endif + + something_was_done = true; + } + + visitor.end_component_handling(); + continue; + } + + // From here on, the CC is a topological disk + if(!remove_self_intersections_with_hole_filling(cc_border_hedges, cc_faces, working_face_range, - tmesh, only_treat_self_intersections_locally, - strong_dihedral_angle, weak_dihedral_angle, - cc_envelope, - vpm, gt)) + tmesh, strong_dihedral_angle, weak_dihedral_angle, + cc_envelope, projector, vpm, gt)) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << " DEBUG: Failed to fill hole\n"; @@ -1912,21 +2271,26 @@ remove_self_intersections_one_step(std::set graph_traits; typedef typename graph_traits::face_descriptor face_descriptor; +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT + IO::write_polygon_mesh("results/input.off", tmesh, parameters::stream_precision(17)); +#endif + // named parameter extraction typedef typename GetVertexPointMap::type VertexPointMap; VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), @@ -1976,12 +2344,15 @@ bool remove_self_intersections(const FaceRange& face_range, GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); bool preserve_genus = choose_parameter(get_parameter(np, internal_np::preserve_genus), true); - const bool only_treat_self_intersections_locally = choose_parameter(get_parameter(np, internal_np::apply_per_connected_component), false); + + // @tmp Squatting that named parameter to signify that treatment should be applied within the CC + // even if there are no self-intersections. For example, two spheres intersecting each other. + const bool treat_all_CCs = choose_parameter(get_parameter(np, internal_np::apply_per_connected_component), true); // When treating intersections locally, we don't want to grow the working range too much as // either the solution is found fast, or it's too difficult and neither local smoothing or local // hole filling are going to provide nice results. - const int default_max_step = only_treat_self_intersections_locally ? 2 : 7; + const int default_max_step = 7; const int max_steps = choose_parameter(get_parameter(np, internal_np::number_of_iterations), default_max_step); // @fixme give it its own named parameter rather than abusing 'with_dihedral_angle'? @@ -1990,7 +2361,8 @@ bool remove_self_intersections(const FaceRange& face_range, // detect_feature_pp NP (unused for now) const double weak_dihedral_angle = 0.; // choose_parameter(get_parameter(np, internal_np::weak_dihedral_angle), 20.); - struct Return_false { + struct Return_false + { bool operator()(std::pair) const { return false; } }; @@ -2005,10 +2377,12 @@ bool remove_self_intersections(const FaceRange& face_range, // use containment check const double containment_epsilon = choose_parameter(get_parameter(np, internal_np::polyhedral_envelope_epsilon), 0.); + internal::Mesh_projection_functor projector(tmesh, vpm); + #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << "DEBUG: Starting remove_self_intersections, is_valid(tmesh)? " << is_valid_polygon_mesh(tmesh) << "\n"; std::cout << "\tpreserve_genus: " << preserve_genus << std::endl; - std::cout << "\tonly_treat_self_intersections_locally: " << only_treat_self_intersections_locally << std::endl; + std::cout << "\ttreat_all_CCs: " << treat_all_CCs << std::endl; std::cout << "\tmax_steps: " << max_steps << std::endl; std::cout << "\tstrong_dihedral_angle: " << strong_dihedral_angle << std::endl; std::cout << "\tweak_dihedral_angle: " << weak_dihedral_angle << std::endl; @@ -2023,7 +2397,7 @@ bool remove_self_intersections(const FaceRange& face_range, Visitor visitor = choose_parameter(get_parameter(np, internal_np::visitor)); visitor.parameters_used(preserve_genus, - only_treat_self_intersections_locally, + treat_all_CCs, max_steps, strong_dihedral_angle, weak_dihedral_angle, @@ -2036,16 +2410,22 @@ bool remove_self_intersections(const FaceRange& face_range, int step = -1; bool all_fixed = true; // indicates if the filling of all created holes went fine bool topology_issue = false; // indicates if some boundary cycles of edges are blocking the fixing - std::set faces_to_remove; + std::set faces_to_treat; std::set working_face_range(face_range.begin(), face_range.end()); visitor.start_main_loop(); while(++step < max_steps) { - if (visitor.stop()) break; +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << " DEBUG: ========== STEP " << step << " / " << max_steps - 1 << " ==========" << std::endl; +#endif + + if(visitor.stop()) + break; + visitor.start_iteration(); - if(faces_to_remove.empty()) // the previous round might have been blocked due to topological constraints + if(faces_to_treat.empty()) // the previous round might have been blocked due to topological constraints { typedef std::pair Face_pair; std::vector self_inter; @@ -2053,32 +2433,33 @@ bool remove_self_intersections(const FaceRange& face_range, // TODO : possible optimization to reduce the range to check with the bbox // of the previous patches or something. self_intersections(working_face_range, tmesh, - CGAL::filter_output_iterator(std::back_inserter(self_inter), out_it_predicates)); + filter_output_iterator(std::back_inserter(self_inter), out_it_predicates), + parameters::vertex_point_map(vpm).geom_traits(gt)); #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG std::cout << " DEBUG: " << self_inter.size() << " intersecting pairs" << std::endl; #endif for(const Face_pair& fp : self_inter) { - faces_to_remove.insert(fp.first); - faces_to_remove.insert(fp.second); + faces_to_treat.insert(fp.first); + faces_to_treat.insert(fp.second); } } - if(faces_to_remove.empty() && all_fixed) + if(faces_to_treat.empty()) { #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG - std::cout << "DEBUG: There are no more faces to remove." << std::endl; + std::cout << "DEBUG: There are no more faces to treat." << std::endl; #endif break; } - visitor.status_update(faces_to_remove); + visitor.status_update(faces_to_treat); std::tie(all_fixed, topology_issue) = internal::remove_self_intersections_one_step( - faces_to_remove, working_face_range, tmesh, - step, preserve_genus, only_treat_self_intersections_locally, - strong_dihedral_angle, weak_dihedral_angle, containment_epsilon, vpm, gt, visitor); + faces_to_treat, working_face_range, tmesh, step, + preserve_genus, treat_all_CCs, strong_dihedral_angle, weak_dihedral_angle, + containment_epsilon, projector, vpm, gt, visitor); #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG if(all_fixed && topology_issue) @@ -2094,14 +2475,21 @@ bool remove_self_intersections(const FaceRange& face_range, std::cout << "solved by unconstrained smoothing: " << internal::self_intersections_solved_by_unconstrained_smoothing << std::endl; std::cout << "solved by constrained hole-filling: " << internal::self_intersections_solved_by_constrained_hole_filling << std::endl; std::cout << "solved by unconstrained hole-filling: " << internal::self_intersections_solved_by_unconstrained_hole_filling << std::endl; - std::cout << "unsolved: " << internal::unsolved_self_intersections << std::endl; + std::cout << "issues during CC treatment: " << internal::unsolved_self_intersections << std::endl; #endif #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT std::ofstream("results/final.off") << std::setprecision(17) << tmesh; #endif - return step < max_steps; + bool self_intersects = does_self_intersect(working_face_range, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)); + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + if(self_intersects) + std::cout << "DEBUG: Failed to solve all self-intersections.\n"; +#endif + + return !self_intersects; } template diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h index 8af96b4cb2e..18035e7f675 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/triangulate_hole.h @@ -136,7 +136,7 @@ namespace Polygon_mesh_processing { CGAL_precondition(face(border_halfedge, pmesh) == boost::graph_traits::null_face()); bool use_cdt = - #ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2 +#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2 false; #else choose_parameter(get_parameter(np, internal_np::use_2d_constrained_delaunay_triangulation), false); @@ -176,22 +176,6 @@ namespace Polygon_mesh_processing { max_squared_distance).first; } - template - void test_in_edges(const PM& pmesh, const VertexRange& patch) - { - for(typename boost::graph_traits::vertex_descriptor v0 : patch) - { - typename boost::graph_traits::in_edge_iterator e, e_end; - for (boost::tie(e, e_end) = in_edges(v0, pmesh); e != e_end; e++) - { - typename boost::graph_traits::halfedge_descriptor he = halfedge(*e, pmesh); - if (is_border(he, pmesh)) { continue; } - - CGAL_assertion(v0 == target(he, pmesh) || v0 == source(he, pmesh)); - } - } - } - /*! \ingroup hole_filling_grp @brief triangulates and refines a hole in a polygon mesh. @@ -285,8 +269,6 @@ namespace Polygon_mesh_processing { triangulate_hole(pmesh, border_halfedge, std::back_inserter(patch), np); face_out = std::copy(patch.begin(), patch.end(), face_out); - test_in_edges(pmesh, vertices(pmesh)); - return refine(pmesh, patch, face_out, vertex_out, np); } @@ -401,16 +383,13 @@ namespace Polygon_mesh_processing { VertexOutputIterator vertex_out, const NamedParameters& np = parameters::default_values()) { + CGAL_precondition(CGAL::is_triangle_mesh(pmesh)); + std::vector::vertex_descriptor> patch; - - CGAL_assertion(CGAL::is_triangle_mesh(pmesh)); - face_out = triangulate_and_refine_hole (pmesh, border_halfedge, face_out, std::back_inserter(patch), np).first; - CGAL_assertion(CGAL::is_triangle_mesh(pmesh)); - - test_in_edges(pmesh, patch); + CGAL_postcondition(CGAL::is_triangle_mesh(pmesh)); bool fair_success = fair(pmesh, patch, np); @@ -536,9 +515,9 @@ bool use_dt3 = #ifndef CGAL_HOLE_FILLING_DO_NOT_USE_CDT2 if (use_cdt) { - struct Always_valid{ - bool operator()(const std::vector&, int,int,int)const - {return true;} + struct Always_valid + { + bool operator()(const std::vector&, int,int,int) const { return true; } }; Always_valid is_valid; @@ -549,8 +528,9 @@ bool use_dt3 = const typename Kernel::FT threshold_distance = choose_parameter( get_parameter(np, internal_np::threshold_distance), typename Kernel::FT(-1)); typename Kernel::FT max_squared_distance = default_squared_distance; - if (threshold_distance >= typename Kernel::FT(0)) + if(threshold_distance >= typename Kernel::FT(0)) max_squared_distance = threshold_distance * threshold_distance; + CGAL_assertion(max_squared_distance >= typename Kernel::FT(0)); if (triangulate_hole_polyline_with_cdt( points, diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h index 8cfafc19b0a..2127f3c8ed3 100644 --- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h +++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h @@ -1529,6 +1529,9 @@ public: /// performs a validity check on a single vertex. bool is_valid(Vertex_index v) const { + if(!has_valid_index(v)) + return false; + Halfedge_index h = vconn_[v].halfedge_; if(h!= null_halfedge() && (!has_valid_index(h) || is_removed(h))) { std::cerr << "Vertex connectivity halfedge error in " << (size_type)v @@ -1540,6 +1543,9 @@ public: /// performs a validity check on a single halfedge. bool is_valid(Halfedge_index h) const { + if(!has_valid_index(h)) + return false; + Face_index f = hconn_[h].face_; Vertex_index v = hconn_[h].vertex_; Halfedge_index hn = hconn_[h].next_halfedge_; @@ -1581,6 +1587,9 @@ public: /// performs a validity check on a single edge. bool is_valid(Edge_index e) const { + if(!has_valid_index(e)) + return false; + Halfedge_index h = halfedge(e); return is_valid(h) && is_valid(opposite(h)); } @@ -1588,6 +1597,9 @@ public: /// performs a validity check on a single face. bool is_valid(Face_index f) const { + if(!has_valid_index(f)) + return false; + Halfedge_index h = fconn_[f].halfedge_; if(!has_valid_index(h) || is_removed(h)) { std::cerr << "Face connectivity halfedge error in " << (size_type)f