diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 95a539040ac..9f9ce51113b 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -86,6 +86,7 @@ create_single_source_cgal_program( "manifoldness_repair_example.cpp" ) create_single_source_cgal_program( "repair_polygon_soup_example.cpp" ) create_single_source_cgal_program( "mesh_smoothing_example.cpp") create_single_source_cgal_program( "locate_example.cpp") +create_single_source_cgal_program( "mesh_self_intersection_removal_example.cpp") if(OpenMesh_FOUND) create_single_source_cgal_program( "compute_normals_example_OM.cpp" ) @@ -118,5 +119,8 @@ find_package(Ceres QUIET) if(TARGET ceres) target_compile_definitions( mesh_smoothing_example PRIVATE CGAL_PMP_USE_CERES_SOLVER ) target_link_libraries( mesh_smoothing_example PRIVATE ceres ) + + target_compile_definitions( mesh_self_intersection_removal_example PRIVATE CGAL_PMP_USE_CERES_SOLVER ) + target_link_libraries( mesh_self_intersection_removal_example PRIVATE ceres ) endif(TARGET ceres) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_self_intersection_removal_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_self_intersection_removal_example.cpp new file mode 100644 index 00000000000..9e83a0ac345 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_self_intersection_removal_example.cpp @@ -0,0 +1,83 @@ +#define CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; + +typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; +typedef typename boost::graph_traits::face_descriptor face_descriptor; + +namespace PMP = CGAL::Polygon_mesh_processing; + +template +void read_mesh(const char* filename, + Mesh& sm) +{ + typedef typename K::Point_3 Point; + + std::ifstream in(filename, std::ios::binary); + if(!in.good()) + { + std::cerr << "Error: can't read file: " << filename << std::endl; + std::exit(1); + } + + std::string fn(filename); + if(fn.substr(fn.find_last_of(".") + 1) == "stl") + { + std::vector points; + std::vector > faces; + CGAL::read_STL(in, points, faces); + + if(!CGAL::Polygon_mesh_processing::orient_polygon_soup(points, faces)) + std::cerr << "W: File does not describe a polygon mesh" << std::endl; + + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, faces, sm); + } + else if(fn.substr(fn.find_last_of(".") + 1) == "off") + { + if(!in || !(in >> sm)) + { + std::cerr << "Error: cannot OFF open mesh\n"; + return; + } + } + else + { + std::cerr << "Unknown file type" << std::endl; + return; + } +} + +int main(int argc, char* argv[]) +{ + const char* filename = (argc > 1) ? argv[1] : "/home/mrouxell/DATA/denser_dinosaur_si.off"; + std::ifstream input(filename); + + Mesh mesh; + read_mesh(filename, mesh); + std::cout << num_vertices(mesh) << " vertices and " << num_faces(mesh) << " faces" << std::endl; + + PMP::remove_degenerate_faces(mesh); + + std::ofstream("in.off") << std::setprecision(17) << mesh; + PMP::remove_self_intersections(mesh); + std::ofstream("out.off") << std::setprecision(17) << mesh; + + std::cout << "Success? " << !PMP::does_self_intersect(mesh) << std::endl; + + return 0; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h index 5bc232135b8..9aa341b4378 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -566,10 +566,8 @@ double approximate_Hausdorff_distance( VertexPointMap vpm_2) { std::vector sample_points; - sample_triangle_mesh( - tm1, - std::back_inserter(sample_points), - np); + sample_triangle_mesh(tm1,std::back_inserter(sample_points),np); + return approximate_Hausdorff_distance(sample_points, tm2, vpm_2); } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h index fb40469e47d..19836618a2b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h @@ -21,7 +21,7 @@ #endif #include -#include +#include #include #include @@ -86,10 +86,18 @@ public: const halfedge_descriptor h = halfedge(e, mesh_); const halfedge_descriptor opp_h = opposite(h, mesh_); - vertex_descriptor v0 = source(h, mesh_); - vertex_descriptor v1 = target(h, mesh_); - vertex_descriptor v2 = target(next(h, mesh_), mesh_); - vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_); + const vertex_descriptor v0 = source(h, mesh_); + const vertex_descriptor v1 = target(h, mesh_); + const vertex_descriptor v2 = target(next(h, mesh_), mesh_); + const vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_); + + // Don't want to flip if the other diagonal already exists + // @todo remeshing can be used to still flip those + std::pair other_hd_already_exists = edge(v2, v3, mesh_); + if(other_hd_already_exists.second) + return false; + + // not local Delaunay := sum of the opposite angles is greater than pi const Point_ref p0 = get(vpmap_, v0); const Point_ref p1 = get(vpmap_, v1); const Point_ref p2 = get(vpmap_, v2); @@ -98,17 +106,7 @@ public: double alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_); double beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_); - // not local Delaunay if the sum of the angles is greater than pi - if(alpha + beta <= CGAL_PI) - return false; - - // Don't want to flip if the other diagonal already exists - // @todo remeshing can be used to still flip those - std::pair other_hd_already_exists = edge(v2, v3, mesh_); - if(other_hd_already_exists.second) - return false; - - return true; + return (alpha + beta > CGAL_PI); } template diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/manifoldness.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/manifoldness.h new file mode 100644 index 00000000000..5ffa00ffeda --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/manifoldness.h @@ -0,0 +1,461 @@ +// Copyright (c) 2015-2019 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Sebastien Loriot, +// Mael Rouxel-Labbé +// +#ifndef CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H +#define CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +/// \ingroup PMP_repairing_grp +/// checks whether a vertex of a polygon mesh is non-manifold. +/// +/// @tparam PolygonMesh a model of `HalfedgeListGraph` +/// +/// @param v a vertex of `pm` +/// @param pm a triangle mesh containing `v` +/// +/// \warning This function has linear runtime with respect to the size of the mesh. +/// +/// \sa `duplicate_non_manifold_vertices()` +/// +/// \return `true` if the vertex is non-manifold, `false` otherwise. +template +bool is_non_manifold_vertex(typename boost::graph_traits::vertex_descriptor v, + const PolygonMesh& pm) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef CGAL::dynamic_halfedge_property_t Halfedge_property_tag; + typedef typename boost::property_map::const_type Visited_halfedge_map; + + // Dynamic pmaps do not have default initialization values (yet) + Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm); + for(halfedge_descriptor h : halfedges(pm)) + put(visited_halfedges, h, false); + + std::size_t incident_null_faces_counter = 0; + for(halfedge_descriptor h : halfedges_around_target(v, pm)) + { + put(visited_halfedges, h, true); + if(CGAL::is_border(h, pm)) + ++incident_null_faces_counter; + } + + if(incident_null_faces_counter > 1) + { + // The vertex is the sole connection between two connected components --> non-manifold + return true; + } + + for(halfedge_descriptor h : halfedges(pm)) + { + if(v == target(h, pm)) + { + // Haven't seen that halfedge yet ==> more than one umbrella incident to 'v' ==> non-manifold + if(!get(visited_halfedges, h)) + return true; + } + } + + return false; +} + +namespace internal { + +template +struct Vertex_collector +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + bool has_old_vertex(const vertex_descriptor v) const { return collections.count(v) != 0; } + void tag_old_vertex(const vertex_descriptor v) + { + CGAL_precondition(!has_old_vertex(v)); + collections[v]; + } + + void collect_vertices(vertex_descriptor v1, vertex_descriptor v2) + { + std::vector& verts = collections[v1]; + if(verts.empty()) + verts.push_back(v1); + verts.push_back(v2); + } + + template + void dump(OutputIterator out) + { + typedef std::pair > Pair_type; + for(const Pair_type& p : collections) + *out++ = p.second; + } + + void dump(Emptyset_iterator) { } + + std::map > collections; +}; + +template +typename boost::graph_traits::vertex_descriptor +create_new_vertex_for_sector(typename boost::graph_traits::halfedge_descriptor sector_begin_h, + typename boost::graph_traits::halfedge_descriptor sector_last_h, + PolygonMesh& pm, + const VPM& vpm, + const ConstraintMap& cmap) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + vertex_descriptor old_vd = target(sector_begin_h, pm); + vertex_descriptor new_vd = add_vertex(pm); + put(vpm, new_vd, get(vpm, old_vd)); + + put(cmap, new_vd, true); + + set_halfedge(new_vd, sector_begin_h, pm); + halfedge_descriptor h = sector_begin_h; + do + { + set_target(h, new_vd, pm); + + if(h == sector_last_h) + break; + else + h = prev(opposite(h, pm), pm); + } + while(h != sector_begin_h); // for safety + CGAL_assertion(h != sector_begin_h); + + return new_vd; +} + +template +std::size_t make_umbrella_manifold(typename boost::graph_traits::halfedge_descriptor h, + PolygonMesh& pm, + internal::Vertex_collector& dmap, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + using parameters::get_parameter; + using parameters::choose_parameter; + + typedef typename GetVertexPointMap::type VertexPointMap; + VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, pm)); + + typedef typename internal_np::Lookup_named_param_def // default (no constraint pmap) + >::type VerticesMap; + VerticesMap cmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), + Constant_property_map(false)); + + std::size_t nb_new_vertices = 0; + + vertex_descriptor old_v = target(h, pm); + put(cmap, old_v, true); // store the duplicates + + // count the number of borders + int border_counter = 0; + halfedge_descriptor ih = h, done = ih, border_h = h; + do + { + if(is_border(ih, pm)) + { + border_h = ih; + ++border_counter; + } + + ih = prev(opposite(ih, pm), pm); + } + while(ih != done); + + bool is_non_manifold_within_umbrella = (border_counter > 1); + if(!is_non_manifold_within_umbrella) + { + const bool first_time_meeting_v = !dmap.has_old_vertex(old_v); + if(first_time_meeting_v) + { + // The star is manifold, so if it is the first time we have met that vertex, + // there is nothing to do, we just keep the same vertex. + set_halfedge(old_v, h, pm); // to ensure halfedge(old_v, pm) stays valid + dmap.tag_old_vertex(old_v); // so that we know we have met old_v already, next time, we'll have to duplicate + } + else + { + // This is not the canonical star associated to 'v'. + // Create a new vertex, and move the whole star to that new vertex + halfedge_descriptor last_h = opposite(next(h, pm), pm); + vertex_descriptor new_v = create_new_vertex_for_sector(h, last_h, pm, vpm, cmap); + dmap.collect_vertices(old_v, new_v); + nb_new_vertices = 1; + } + } + // if there is more than one sector, look at each sector and split them away from the main one + else + { + // the first manifold sector, described by two halfedges + halfedge_descriptor sector_start_h = border_h; + CGAL_assertion(is_border(border_h, pm)); + + bool should_stop = false; + bool is_main_sector = true; + do + { + CGAL_assertion(is_border(sector_start_h, pm)); + + // collect the sector and split it away if it must be + halfedge_descriptor sector_last_h = sector_start_h; + do + { + halfedge_descriptor next_h = prev(opposite(sector_last_h, pm), pm); + + if(is_border(next_h, pm)) + break; + + sector_last_h = next_h; + } + while(sector_last_h != sector_start_h); + CGAL_assertion(!is_border(sector_last_h, pm)); + CGAL_assertion(sector_last_h != sector_start_h); + + halfedge_descriptor next_start_h = prev(opposite(sector_last_h, pm), pm); + + // there are multiple CCs incident to this particular vertex, and we should create a new vertex + // if it's not the first umbrella around 'old_v' or not the first sector, but only not if it's + // both the first umbrella and first sector. + bool must_create_new_vertex = (!is_main_sector || dmap.has_old_vertex(old_v)); + + // In any case, we must set up the next pointer correctly + set_next(sector_start_h, opposite(sector_last_h, pm), pm); + + if(must_create_new_vertex) + { + vertex_descriptor new_v = create_new_vertex_for_sector(sector_start_h, sector_last_h, pm, vpm, cmap); + dmap.collect_vertices(old_v, new_v); + ++nb_new_vertices; + } + else + { + // We are in the first sector and first star, ensure that halfedge(old_v, pm) stays valid + set_halfedge(old_v, sector_start_h, pm); + } + + is_main_sector = false; + sector_start_h = next_start_h; + should_stop = (sector_start_h == border_h); + } + while(!should_stop); + } + + return nb_new_vertices; +} + +} // end namespace internal + +/// \ingroup PMP_repairing_grp +/// collects the non-manifold vertices (if any) present in the mesh. A non-manifold vertex `v` is returned +/// via one incident halfedge `h` such that `target(h, pm) = v` for all the umbrellas that `v` apppears in +/// (an umbrella being the set of faces incident to all the halfedges reachable by walking around `v` +/// using `hnext = prev(opposite(h, pm), pm)`, starting from `h`). +/// +/// @tparam PolygonMesh a model of `HalfedgeListGraph` +/// @tparam OutputIterator a model of `OutputIterator` holding objects of type +/// `boost::graph_traits::%halfedge_descriptor` +/// +/// @param pm a triangle mesh +/// @param out the output iterator that collects halfedges incident to `v` +/// +/// \sa `is_non_manifold_vertex()` +/// \sa `duplicate_non_manifold_vertices()` +/// +/// \return the output iterator. +template +OutputIterator non_manifold_vertices(const PolygonMesh& pm, + OutputIterator out) +{ + // 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. + + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef CGAL::dynamic_vertex_property_t Vertex_bool_tag; + typedef typename boost::property_map::const_type Known_manifold_vertex_map; + typedef CGAL::dynamic_vertex_property_t Vertex_halfedge_tag; + typedef typename boost::property_map::const_type Visited_vertex_map; + typedef CGAL::dynamic_halfedge_property_t Halfedge_property_tag; + typedef typename boost::property_map::const_type Visited_halfedge_map; + + Known_manifold_vertex_map known_nm_vertices = get(Vertex_bool_tag(), pm); + Visited_vertex_map visited_vertices = get(Vertex_halfedge_tag(), pm); + Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm); + + halfedge_descriptor null_h = boost::graph_traits::null_halfedge(); + + // Dynamic pmaps do not have default initialization values (yet) + for(vertex_descriptor v : vertices(pm)) + { + put(known_nm_vertices, v, false); + put(visited_vertices, v, null_h); + } + for(halfedge_descriptor h : halfedges(pm)) + put(visited_halfedges, h, false); + + for(halfedge_descriptor h : halfedges(pm)) + { + // If 'h' is not visited yet, we walk around the target of 'h' and mark these + // halfedges as visited. Thus, if we are here and the target is already marked as visited, + // it means that the vertex is non manifold. + if(!get(visited_halfedges, h)) + { + put(visited_halfedges, h, true); + bool is_non_manifold = false; + + vertex_descriptor v = target(h, pm); + + if(get(visited_vertices, v) != null_h) // already seen this vertex, but not from this star + { + is_non_manifold = true; + + // if this is the second time we visit that vertex and the first star was manifold, we have + // never reported the first star, but we must now + if(!get(known_nm_vertices, v)) + *out++ = get(visited_vertices, v); // that's a halfedge of the first star we've seen 'v' in + } + else + { + // first time we meet this vertex, just mark it so, and keep the halfedge we found the vertex with in memory + put(visited_vertices, v, h); + } + + // While walking the star of this halfedge, if we meet a border halfedge more than once, + // it means the mesh is pinched and we are also in the case of a non-manifold situation + halfedge_descriptor ih = h, done = ih; + int border_counter = 0; + do + { + put(visited_halfedges, ih, true); + if(is_border(ih, pm)) + ++border_counter; + + ih = prev(opposite(ih, pm), pm); + } + while(ih != done); + + if(border_counter > 1) + is_non_manifold = true; + + if(is_non_manifold) + { + *out++ = h; + put(known_nm_vertices, v, true); + } + } + } + + return out; +} + +/// \ingroup PMP_repairing_grp +/// duplicates all the non-manifold vertices of the input mesh. +/// +/// @tparam PolygonMesh a model of `HalfedgeListGraph` and `MutableHalfedgeGraph` +/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" +/// +/// @param pm the surface mesh to be repaired +/// @param np optional \ref pmp_namedparameters "Named Parameters" described below +/// +/// \cgalNamedParamsBegin +/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. +/// The type of this map is model of `ReadWritePropertyMap`. +/// If this parameter is omitted, an internal property map for +/// `CGAL::vertex_point_t` should be available in `PolygonMesh` +/// \cgalParamEnd +/// \cgalParamBegin{vertex_is_constrained_map} a writable property map with `vertex_descriptor` +/// as key and `bool` as `value_type`. `put(pmap, v, true)` will be called for each duplicated +/// vertices, as well as the original non-manifold vertex in the input mesh. +/// \cgalParamEnd +/// \cgalParamBegin{output_iterator} a model of `OutputIterator` with value type +/// `std::vector`. The first vertex of each vector is a non-manifold vertex +/// of the input mesh, followed by the new vertices that were created to fix this precise +/// non-manifold configuration. +/// \cgalParamEnd +/// \cgalNamedParamsEnd +/// +/// \return the number of vertices created. +template +std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm, + const NamedParameters& np) +{ + using parameters::get_parameter; + using parameters::choose_parameter; + + typedef boost::graph_traits GT; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + + typedef typename internal_np::Lookup_named_param_def::type Output_iterator; + + Output_iterator out = choose_parameter(get_parameter(np, internal_np::output_iterator), + Emptyset_iterator()); + + std::vector non_manifold_cones; + non_manifold_vertices(pm, std::back_inserter(non_manifold_cones)); + + internal::Vertex_collector dmap; + std::size_t nb_new_vertices = 0; + if(!non_manifold_cones.empty()) + { + for(halfedge_descriptor h : non_manifold_cones) + nb_new_vertices += internal::make_umbrella_manifold(h, pm, dmap, np); + + dmap.dump(out); + } + + return nb_new_vertices; +} + +template +std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm) +{ + return duplicate_non_manifold_vertices(pm, parameters::all_default()); +} + +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#endif // CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_degeneracies.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_degeneracies.h new file mode 100644 index 00000000000..a61062ddf25 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_degeneracies.h @@ -0,0 +1,1814 @@ +// Copyright (c) 2015-2019 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Sebastien Loriot, +// Mael Rouxel-Labbé +// +#ifndef CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H +#define CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H + +#include + +#include + +#include +#include + +#include + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG +#include +#include +#endif + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { +namespace debug { + +template +std::ostream& +dump_edge_neighborhood(typename boost::graph_traits::edge_descriptor ed, + TriangleMesh& tmesh, + const VertexPointMap& vpmap, + std::ostream& out) +{ + typedef boost::graph_traits GT; + typedef typename GT::vertex_descriptor vertex_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::face_descriptor face_descriptor; + + halfedge_descriptor h = halfedge(ed, tmesh); + + std::map vertices; + std::set faces; + int vindex = 0; + + for(halfedge_descriptor hd : halfedges_around_target(h, tmesh)) + { + if(vertices.insert(std::make_pair(source(hd, tmesh), vindex)).second) + ++vindex; + + if(!is_border(hd, tmesh)) + faces.insert(face(hd, tmesh)); + } + + h = opposite(h, tmesh); + for(halfedge_descriptor hd : halfedges_around_target(h, tmesh)) + { + if(vertices.insert(std::make_pair(source(hd, tmesh), vindex)).second) + ++vindex; + + if(!is_border(hd, tmesh)) + faces.insert(face(hd, tmesh)); + } + + std::vector ordered_vertices(vertices.size()); + typedef std::pair Pair_type; + for(const Pair_type& p : vertices) + ordered_vertices[p.second] = p.first; + + out << "OFF\n" << ordered_vertices.size() << " " << faces.size() << " 0\n"; + for(vertex_descriptor vd : ordered_vertices) + out << get(vpmap, vd) << "\n"; + for(face_descriptor fd : faces) + { + out << "3"; + h = halfedge(fd, tmesh); + for(halfedge_descriptor hd : halfedges_around_face(h, tmesh)) + out << " " << vertices[target(hd, tmesh)]; + out << "\n"; + } + + return out; +} + +template +void dump_cc_faces(const FaceRange& cc_faces, const TriangleMesh& tm, std::ostream& output) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename boost::property_map::const_type Vpm; + typedef typename boost::property_traits::value_type Point_3; + + Vpm vpm = get(boost::vertex_point, tm); + + int id = 0; + std::map vids; + + for(face_descriptor f : cc_faces) + { + if(vids.insert(std::make_pair(target(halfedge(f, tm), tm), id)).second) ++id; + if(vids.insert(std::make_pair(target(next(halfedge(f, tm), tm), tm), id)).second) ++id; + if(vids.insert(std::make_pair(target(next(next(halfedge(f, tm), tm), tm), tm), id)).second) ++id; + } + + std::vector points(vids.size()); + typedef std::pair Pair_type; + for(const Pair_type& p : vids) + points[p.second] = get(vpm, p.first); + + output << std::setprecision(17); + output << "OFF\n" << vids.size() << " " << cc_faces.size() << " 0\n"; + for(Point_3 p : points) + output << p << "\n"; + for(face_descriptor f : cc_faces) + { + output << "3 " + << vids[target(halfedge(f, tm), tm)] << " " + << vids[target(next(halfedge(f, tm), tm), tm)] << " " + << vids[target(next(next(halfedge(f, tm), tm), tm), tm)] << "\n"; + } +} + +} // namespace debug + +namespace internal { + +template +struct Less_vertex_point +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + Less_vertex_point(const Traits& traits, const VertexPointMap& vpmap) + : m_traits(traits), m_vpmap(vpmap) + {} + + bool operator()(vertex_descriptor v1, vertex_descriptor v2) const { + return m_traits.less_xyz_3_object()(get(m_vpmap, v1), get(m_vpmap, v2)); + } + + const Traits& m_traits; + const VertexPointMap& m_vpmap; +}; + +} // end namespace internal + +// this function removes a border edge even if it does not satisfy the link condition. +// null_vertex() is returned if the removal changes the topology of the input +template +typename boost::graph_traits::vertex_descriptor +remove_a_border_edge(typename boost::graph_traits::edge_descriptor ed, + TriangleMesh& tm, + EdgeSet& edge_set, + FaceSet& face_set) +{ + typedef boost::graph_traits GT; + typedef typename GT::vertex_descriptor vertex_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::edge_descriptor edge_descriptor; + typedef typename GT::face_descriptor face_descriptor; + + halfedge_descriptor h = halfedge(ed, tm); + + if(is_border(h, tm)) + h = opposite(h, tm); + + halfedge_descriptor opp_h = opposite(h, tm); + CGAL_assertion(is_border(opp_h, tm)); + CGAL_assertion(!is_border(h, tm)); + + CGAL_assertion(next(next(opp_h, tm), tm) != opp_h); // not working for a hole made of 2 edges + CGAL_assertion(next(next(next(opp_h, tm), tm), tm) != opp_h); // not working for a hole make of 3 edges + + if(CGAL::Euler::does_satisfy_link_condition(edge(h, tm), tm)) + { + edge_set.erase(ed); + halfedge_descriptor h = halfedge(ed, tm); + if(is_border(h, tm)) + h = opposite(h, tm); + + edge_set.erase(edge(prev(h, tm), tm)); + face_set.erase(face(h, tm)); + + return CGAL::Euler::collapse_edge(ed, tm); + } + + // collect edges that have one vertex in the link of + // the vertices of h and one of the vertex of h as other vertex + std::set common_incident_edges; + for(halfedge_descriptor hos : halfedges_around_source(h, tm)) + { + for(halfedge_descriptor hot : halfedges_around_target(h, tm)) + { + if(target(hos, tm) == source(hot, tm)) + { + common_incident_edges.insert(edge(hot, tm)); + common_incident_edges.insert(edge(hos, tm)); + } + } + } + + CGAL_assertion(common_incident_edges.size() >= 2); + + // in the following loop, we visit define a connected component of + // faces bounded by edges in common_incident_edges and h. We look + // for the maximal one. This set of faces is the one that will + // disappear while collapsing ed + std::set marked_faces; + + std::vector queue; + queue.push_back(opposite(next(h, tm), tm)); + queue.push_back(opposite(prev(h, tm), tm)); + marked_faces.insert(face(h, tm)); + + do + { + std::vector boundary; + while(!queue.empty()) + { + halfedge_descriptor back=queue.back(); + queue.pop_back(); + face_descriptor fback=face(back, tm); + if(common_incident_edges.count(edge(back, tm))) + { + boundary.push_back(back); + continue; + } + + if(fback==GT::null_face() || !marked_faces.insert(fback).second) + continue; + + queue.push_back(opposite(next(back, tm), tm)); + if(is_border(queue.back(), tm)) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Boundary reached during exploration, the region to be removed is not a topological disk, not handled for now.\n"; +#endif + return GT::null_vertex(); + } + + queue.push_back(opposite(prev(back, tm), tm)); + if(is_border(queue.back(), tm)) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Boundary reached during exploration, the region to be removed is not a topological disk, not handled for now.\n"; +#endif + return GT::null_vertex(); + } + } + + CGAL_assertion(boundary.size() == 2); + common_incident_edges.erase(edge(boundary[0], tm)); + common_incident_edges.erase(edge(boundary[1], tm)); + if(!is_border(boundary[0], tm) || common_incident_edges.empty()) + queue.push_back(boundary[0]); + if(!is_border(boundary[1], tm) || common_incident_edges.empty()) + queue.push_back(boundary[1]); + } + while(!common_incident_edges.empty()); + + // hk1 and hk2 are bounding the region that will be removed. + // The edge of hk2 will be removed and hk2 will be replaced + // by the opposite edge of hk1 + halfedge_descriptor hk1 = queue.front(); + halfedge_descriptor hk2 = queue.back(); + if(target(hk1, tm)!=source(hk2, tm)) + std::swap(hk1, hk2); + + CGAL_assertion(target(hk1, tm) == source(hk2, tm)); + CGAL_assertion(source(hk1, tm) == source(h, tm)); + CGAL_assertion(target(hk2, tm) == target(h, tm)); + + CGAL_assertion(is_valid_polygon_mesh(tm)); + if(!is_selection_a_topological_disk(marked_faces, tm)) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "The region to be removed is not a topological disk, not handled for now.\n"; +#endif + return GT::null_vertex(); + } + + if(is_border(hk1, tm) && is_border(hk2, tm)) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "The region to be removed is an isolated region, not handled for now.\n"; +#endif + return GT::null_vertex(); + } + + // collect vertices and edges to remove and do remove faces + std::set edges_to_remove; + std::set vertices_to_remove; + for(face_descriptor fd : marked_faces) + { + halfedge_descriptor hd = halfedge(fd, tm); + for(int i=0; i<3; ++i) + { + edges_to_remove.insert(edge(hd, tm)); + vertices_to_remove.insert(target(hd, tm)); + hd = next(hd, tm); + } + } + + vertex_descriptor vkept = source(hk1, tm); + + //back-up next, prev halfedge pointers to be restored after removal + halfedge_descriptor hp = prev(opp_h, tm); + halfedge_descriptor hn = next(opp_h, tm); + halfedge_descriptor hk1_opp_next = next(hk2, tm); + halfedge_descriptor hk1_opp_prev = prev(hk2, tm); + face_descriptor hk1_opp_face = face(hk2, tm); + + // we will remove the target of hk2, update vertex pointers + for(halfedge_descriptor hot : halfedges_around_target(hk2, tm)) + set_target(hot, vkept, tm); + + // update halfedge pointers since hk2 will be removed + set_halfedge(vkept, opposite(hk1, tm), tm); + set_halfedge(target(hk1, tm), hk1, tm); + + // do not remove hk1 and its vertices + vertices_to_remove.erase(vkept); + vertices_to_remove.erase(target(hk1, tm)); + edges_to_remove.erase(edge(hk1, tm)); + + bool hk2_equals_hp = (hk2 == hp); + CGAL_assertion(is_border(hk2, tm) == hk2_equals_hp); + + /* + - case hk2!=hp + + /\ / + hk1/ \hk2 / + / \ / + ____/______\/____ + hn h_opp hp + + - case hk2==hp + + /\ + hk1/ \hk2 == hp + / \ + ____/______\ + hn h_opp + */ + + // remove vertices + for(vertex_descriptor vd : vertices_to_remove) + remove_vertex(vd, tm); + + // remove edges + for(edge_descriptor ed : edges_to_remove) + { + edge_set.erase(ed); + remove_edge(ed, tm); + } + + // remove faces + for(face_descriptor fd : marked_faces) + { + face_set.erase(fd); + remove_face(fd, tm); + } + + // now update pointers + set_face(opposite(hk1, tm), hk1_opp_face, tm); + if(!hk2_equals_hp) + { + set_next(hp, hn, tm); + set_next(opposite(hk1, tm), hk1_opp_next, tm); + set_next(hk1_opp_prev, opposite(hk1, tm), tm); + set_halfedge(hk1_opp_face, opposite(hk1, tm), tm); + } + else + { + set_next(hk1_opp_prev, opposite(hk1, tm), tm); + set_next(opposite(hk1, tm), hn, tm); + } + + CGAL_assertion(is_valid_polygon_mesh(tm)); + return vkept; +} + +template +typename boost::graph_traits::vertex_descriptor +remove_a_border_edge(typename boost::graph_traits::edge_descriptor ed, + TriangleMesh& tm) +{ + std::set::edge_descriptor> edge_set; + std::set::face_descriptor> face_set; + + return remove_a_border_edge(ed, tm, edge_set, face_set); +} + +template +bool remove_degenerate_edges(const EdgeRange& edge_range, + TriangleMesh& tmesh, + FaceSet& face_set, + const NamedParameters& np) +{ + CGAL_assertion(CGAL::is_triangle_mesh(tmesh)); + CGAL_assertion(CGAL::is_valid_polygon_mesh(tmesh)); + + using parameters::get_parameter; + using parameters::choose_parameter; + + typedef TriangleMesh TM; + typedef typename boost::graph_traits GT; + typedef typename GT::vertex_descriptor vertex_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::edge_descriptor edge_descriptor; + typedef typename GT::face_descriptor face_descriptor; + + typedef typename GetVertexPointMap::type VertexPointMap; + VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, tmesh)); + + typedef typename GetGeomTraits::type Traits; + + std::size_t nb_deg_faces = 0; + bool all_removed = false; + bool some_removed = true; + bool preserve_genus = choose_parameter(get_parameter(np, internal_np::preserve_genus), true); + + // collect edges of length 0 + while(some_removed && !all_removed) + { + some_removed = false; + all_removed = true; + std::set degenerate_edges_to_remove; + degenerate_edges(edge_range, tmesh, std::inserter(degenerate_edges_to_remove, + degenerate_edges_to_remove.end())); + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Found " << degenerate_edges_to_remove.size() << " null edges.\n"; +#endif + + // first try to remove all collapsable edges + typename std::set::iterator it = degenerate_edges_to_remove.begin(); + while(it != degenerate_edges_to_remove.end()) + { + edge_descriptor e = *it; + if(CGAL::Euler::does_satisfy_link_condition(e, tmesh)) + { + halfedge_descriptor h = halfedge(e, tmesh); + degenerate_edges_to_remove.erase(it); + + // remove edges that could also be set for removal + if(face(h, tmesh) != GT::null_face()) + { + ++nb_deg_faces; + degenerate_edges_to_remove.erase(edge(prev(h, tmesh), tmesh)); + face_set.erase(face(h, tmesh)); + } + + if(face(opposite(h, tmesh), tmesh) != GT::null_face()) + { + ++nb_deg_faces; + degenerate_edges_to_remove.erase(edge(prev(opposite(h, tmesh), tmesh), tmesh)); + face_set.erase(face(opposite(h, tmesh), tmesh)); + } + + //now remove the edge + CGAL::Euler::collapse_edge(e, tmesh); + + // some_removed is not updated on purpose because if nothing + // happens below then nothing can be done + it = degenerate_edges_to_remove.begin(); + } + else // link condition not satisfied + { + ++it; + } + } + + CGAL_assertion(is_valid_polygon_mesh(tmesh)); +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Remaining " << degenerate_edges_to_remove.size() << " null edges to be handled.\n"; +#endif + + while(!degenerate_edges_to_remove.empty()) + { + edge_descriptor ed = *degenerate_edges_to_remove.begin(); + degenerate_edges_to_remove.erase(degenerate_edges_to_remove.begin()); + halfedge_descriptor h = halfedge(ed, tmesh); + + if(CGAL::Euler::does_satisfy_link_condition(ed, tmesh)) + { + // remove edges that could also be set for removal + if(face(h, tmesh) != GT::null_face()) + { + ++nb_deg_faces; + degenerate_edges_to_remove.erase(edge(prev(h, tmesh), tmesh)); + face_set.erase(face(h, tmesh)); + } + + if(face(opposite(h, tmesh), tmesh)!=GT::null_face()) + { + ++nb_deg_faces; + degenerate_edges_to_remove.erase(edge(prev(opposite(h, tmesh), tmesh), tmesh)); + face_set.erase(face(opposite(h, tmesh), tmesh)); + } + + //now remove the edge + CGAL::Euler::collapse_edge(ed, tmesh); + some_removed = true; + } + else // link condition not satisfied + { + // handle the case when the edge is incident to a triangle hole + // we first fill the hole and try again + if(is_border(ed, tmesh)) + { + halfedge_descriptor hd = halfedge(ed, tmesh); + if(!is_border(hd, tmesh)) + hd = opposite(hd, tmesh); + + if(is_triangle(hd, tmesh)) + { + if(!preserve_genus) + { + Euler::fill_hole(hd, tmesh); + degenerate_edges_to_remove.insert(ed); // reinsert the edge for future processing + } + else + { + all_removed=false; + } + + continue; + } + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Calling remove_a_border_edge\n"; +#endif + + vertex_descriptor vd = remove_a_border_edge(ed, tmesh, degenerate_edges_to_remove, face_set); + if(vd == GT::null_vertex()) + { + // TODO: if some border edges are later removed, the edge might be processable later + // for example if it belongs to boundary cycle of edges where the number of non-degenerate + // edges is 2. That's what happen with fused_vertices.off in the testsuite where the edges + // are not processed the same way with Polyhedron and Surface_mesh. In the case of Polyhedron + // more degenerate edges could be removed. + all_removed = false; + } + else + some_removed = true; + + continue; + } + else + { + halfedge_descriptor hd = halfedge(ed, tmesh); + // if both vertices are boundary vertices we can't do anything + bool impossible = false; + for(halfedge_descriptor h : halfedges_around_target(hd, tmesh)) + { + if(is_border(h, tmesh)) + { + impossible = true; + break; + } + } + + if(impossible) + { + impossible = false; + for(halfedge_descriptor h : halfedges_around_source(hd, tmesh)) + { + if(is_border(h, tmesh)) + { + impossible = true; + break; + } + } + + if(impossible) + { + all_removed = false; + continue; + } + } + } + + // When the edge does not satisfy the link condition, it means that it cannot be + // collapsed as is. In the following if there is a topological issue + // with contracting the edge (component or geometric feature that disappears), + // nothing is done. + // We start by marking the faces that are incident to an edge endpoint. + // If the set of marked faces is a topologically disk, then we simply remove all the simplicies + // inside the disk and star the hole with the edge vertex kept. + // If the set of marked faces is not a topological disk, it has some non-manifold vertices + // on its boundary. We need to mark additional faces to make it a topological disk. + // We can then apply the star hole procedure. + // Right now we additionally mark the smallest connected components of non-marked faces + // (using the number of faces) + + //backup central point + typename Traits::Point_3 pt = get(vpmap, source(ed, tmesh)); + + // mark faces of the link of each endpoints of the edge which collapse is not topologically valid + std::set marked_faces; + + // first endpoint + for(halfedge_descriptor hd : CGAL::halfedges_around_target(halfedge(ed, tmesh), tmesh)) + if(!is_border(hd, tmesh)) + marked_faces.insert(face(hd, tmesh)); + + // second endpoint + for(halfedge_descriptor hd : CGAL::halfedges_around_target(opposite(halfedge(ed, tmesh), tmesh), tmesh)) + if(!is_border(hd, tmesh)) + marked_faces.insert(face(hd, tmesh)); + + // extract the halfedges on the boundary of the marked region + std::vector border; + for(face_descriptor fd : marked_faces) + { + for(halfedge_descriptor hd : CGAL::halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + halfedge_descriptor hd_opp = opposite(hd, tmesh); + if(is_border(hd_opp, tmesh) || marked_faces.count(face(hd_opp, tmesh)) == 0) + { + border.push_back(hd); + } + } + } + + if(border.empty()) + { + // a whole connected component (without boundary) got selected and will disappear (not handled for now) +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Trying to remove a whole connected component, not handled yet\n"; +#endif + all_removed = false; + continue; + } + + // define cc of border halfedges: two halfedges are in the same cc + // if they are on the border of the cc of non-marked faces. + typedef CGAL::Union_find UF_ds; + UF_ds uf; + std::map handles; + + // one cc per border halfedge + for(halfedge_descriptor hd : border) + handles.insert(std::make_pair(hd, uf.make_set(hd))); + + // join cc's + for(halfedge_descriptor hd : border) + { + CGAL_assertion(marked_faces.count(face(hd, tmesh)) > 0); + CGAL_assertion(marked_faces.count(face(opposite(hd, tmesh), tmesh)) == 0); + halfedge_descriptor candidate = hd; + + do + { + candidate = prev(opposite(candidate, tmesh), tmesh); + } + while(!marked_faces.count(face(opposite(candidate, tmesh), tmesh))); + + uf.unify_sets(handles[hd], handles[opposite(candidate, tmesh)]); + } + + std::size_t nb_cc = uf.number_of_sets(); + if(nb_cc != 1) + { + // if more than one connected component is found then the patch + // made of marked faces contains "non-manifold" vertices. + // The smallest components need to be marked so that the patch + // made of marked faces is a topological disk + + // we will explore in parallel the connected components and will stop + // when all but one connected component have been entirely explored. + // We add one face at a time for each cc in order to not explore a + // potentially very large cc. + std::vector > stacks_per_cc(nb_cc); + std::vector > faces_per_cc(nb_cc); + std::vector exploration_finished(nb_cc, false); + + // init the stacks of halfedges using the cc of the boundary + std::size_t index = 0; + std::map ccs; + + typedef std::pair Pair_type; + for(const Pair_type& p : handles) + { + halfedge_descriptor opp_hedge = opposite(p.first, tmesh); + if(is_border(opp_hedge, tmesh)) + continue; // nothing to do on the boundary + + typedef typename std::map::iterator Map_it; + std::pair insert_res = ccs.insert(std::make_pair(*uf.find(p.second), index)); + + if(insert_res.second) + ++index; + + stacks_per_cc[insert_res.first->second].push_back(prev(opp_hedge, tmesh)); + stacks_per_cc[insert_res.first->second].push_back(next(opp_hedge, tmesh)); + faces_per_cc[insert_res.first->second].insert(face(opp_hedge, tmesh)); + } + + if(index != nb_cc) + { + // most probably, one cc is a cycle of border edges +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Trying to remove a component with a cycle of halfedges (nested hole or whole component), not handled yet.\n"; +#endif + all_removed = false; + continue; + } + + std::size_t nb_ccs_to_be_explored = nb_cc; + index = 0; + //explore the cc's + do + { + // try to extract one more face for a given cc + do + { + CGAL_assertion(!exploration_finished[index]); + CGAL_assertion(!stacks_per_cc.empty()); + CGAL_assertion(!stacks_per_cc[index].empty()); + + halfedge_descriptor hd = stacks_per_cc[index].back(); + stacks_per_cc[index].pop_back(); + hd = opposite(hd, tmesh); + + if(!is_border(hd, tmesh) && !marked_faces.count(face(hd, tmesh))) + { + if(faces_per_cc[index].insert(face(hd, tmesh)).second) + { + stacks_per_cc[index].push_back(next(hd, tmesh)); + stacks_per_cc[index].push_back(prev(hd, tmesh)); + break; + } + } + + if(stacks_per_cc[index].empty()) + break; + } + while(true); + + // the exploration of a cc is finished when its stack is empty + exploration_finished[index]=stacks_per_cc[index].empty(); + if(exploration_finished[index]) + --nb_ccs_to_be_explored; + + if(nb_ccs_to_be_explored==1) + break; + + while(exploration_finished[(++index)%nb_cc]); + + index = index%nb_cc; + } + while(true); + + /// \todo use the area criteria? this means maybe continue exploration of larger cc + // mark faces of completetly explored cc + for(index=0; index vertices_to_keep; + std::set halfedges_to_keep; + for(halfedge_descriptor hd : border) + { + if(!marked_faces.count(face(opposite(hd, tmesh), tmesh))) + { + halfedges_to_keep.insert(hd); + vertices_to_keep.insert(target(hd, tmesh)); + } + } + + // backup next,prev relationships to set after patch removal + std::vector > next_prev_halfedge_pairs; + halfedge_descriptor first_border_hd = *(halfedges_to_keep.begin()); + halfedge_descriptor current_border_hd = first_border_hd; + do + { + halfedge_descriptor prev_border_hd = current_border_hd; + current_border_hd = next(current_border_hd, tmesh); + + while(marked_faces.count(face(opposite(current_border_hd, tmesh), tmesh))) + current_border_hd=next(opposite(current_border_hd, tmesh), tmesh); + + next_prev_halfedge_pairs.push_back(std::make_pair(prev_border_hd, current_border_hd)); + } + while(current_border_hd!=first_border_hd); + + // collect vertices and edges to remove and do remove faces + std::set edges_to_remove; + std::set vertices_to_remove; + for(face_descriptor fd : marked_faces) + { + halfedge_descriptor hd=halfedge(fd, tmesh); + for(int i=0; i<3; ++i) + { + if(!halfedges_to_keep.count(hd)) + edges_to_remove.insert(edge(hd, tmesh)); + + if(!vertices_to_keep.count(target(hd, tmesh))) + vertices_to_remove.insert(target(hd, tmesh)); + + hd = next(hd, tmesh); + } + + remove_face(fd, tmesh); + face_set.erase(fd); + } + + // remove vertices + for(vertex_descriptor vd : vertices_to_remove) + remove_vertex(vd, tmesh); + + // remove edges + for(edge_descriptor ed : edges_to_remove) + { + degenerate_edges_to_remove.erase(ed); + remove_edge(ed, tmesh); + } + + // add a new face, set all border edges pointing to it + // and update halfedge vertex of patch boundary vertices + face_descriptor new_face = add_face(tmesh); + typedef std::pair Pair_type; + for(const Pair_type& p : next_prev_halfedge_pairs) + { + set_face(p.first, new_face, tmesh); + set_next(p.first, p.second, tmesh); + set_halfedge(target(p.first, tmesh), p.first, tmesh); + } + + set_halfedge(new_face, first_border_hd, tmesh); + + // triangulate the new face and update the coordinate of the central vertex + halfedge_descriptor new_hd=Euler::add_center_vertex(first_border_hd, tmesh); + put(vpmap, target(new_hd, tmesh), pt); + + for(halfedge_descriptor hd : halfedges_around_target(new_hd, tmesh)) + { + if(is_degenerate_edge(edge(hd, tmesh), tmesh, np)) + degenerate_edges_to_remove.insert(edge(hd, tmesh)); + + if(face(hd, tmesh) != GT::null_face() && is_degenerate_triangle_face(face(hd, tmesh), tmesh)) + face_set.insert(face(hd, tmesh)); + } + + CGAL_assertion(is_valid_polygon_mesh(tmesh)); + } + } + } + + return all_removed; +} + +template +bool remove_degenerate_edges(const EdgeRange& edge_range, + TriangleMesh& tmesh, + const CGAL_PMP_NP_CLASS& np) +{ + std::set::face_descriptor> face_set; + return remove_degenerate_edges(edge_range, tmesh, face_set, np); +} + +template +bool remove_degenerate_edges(TriangleMesh& tmesh, + const CGAL_PMP_NP_CLASS& np) +{ + std::set::face_descriptor> face_set; + return remove_degenerate_edges(edges(tmesh), tmesh, face_set, np); +} + +template +bool remove_degenerate_edges(const EdgeRange& edge_range, + TriangleMesh& tmesh) +{ + std::set::face_descriptor> face_set; + return remove_degenerate_edges(edge_range, tmesh, face_set, parameters::all_default()); +} + +template +bool remove_degenerate_edges(TriangleMesh& tmesh) +{ + std::set::face_descriptor> face_set; + return remove_degenerate_edges(edges(tmesh), tmesh, face_set, parameters::all_default()); +} + +// \ingroup PMP_repairing_grp +// removes the degenerate faces from a triangulated surface mesh. +// A face is considered degenerate if two of its vertices share the same location, +// or more generally if all its vertices are collinear. +// +// @pre `CGAL::is_triangle_mesh(tmesh)` +// +// @tparam TriangleMesh a model of `FaceListGraph` and `MutableFaceGraph` +// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" +// +// @param tmesh the triangulated surface mesh to be repaired +// @param np optional \ref pmp_namedparameters "Named Parameters" described below +// +// \cgalNamedParamsBegin +// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. +// The type of this map is model of `ReadWritePropertyMap`. +// If this parameter is omitted, an internal property map for +// `CGAL::vertex_point_t` must be available in `TriangleMesh` +// \cgalParamEnd +// \cgalParamBegin{geom_traits} a geometric traits class instance. +// The traits class must provide the nested type `Point_3`, +// and the nested functors: +// - `Compare_distance_3` to compute the distance between 2 points +// - `Collinear_3` to check whether 3 points are collinear +// - `Less_xyz_3` to compare lexicographically two points +// - `Equal_3` to check whether 2 points are identical. +// For each functor Foo, a function `Foo foo_object()` must be provided. +// \cgalParamEnd +// \cgalNamedParamsEnd +// +// \todo the function might not be able to remove all degenerate faces. +// We should probably do something with the return type. +// +// \return `true` if all degenerate faces were successfully removed, and `false` otherwise. +template +bool remove_degenerate_faces(const FaceRange& face_range, + TriangleMesh& tmesh, + const NamedParameters& np) +{ + CGAL_assertion(CGAL::is_triangle_mesh(tmesh)); + + using parameters::get_parameter; + using parameters::choose_parameter; + + typedef TriangleMesh TM; + typedef typename boost::graph_traits GT; + typedef typename GT::vertex_descriptor vertex_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::edge_descriptor edge_descriptor; + typedef typename GT::face_descriptor face_descriptor; + + typedef typename GetVertexPointMap::type VertexPointMap; + VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, tmesh)); + typedef typename GetGeomTraits::type Traits; + Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits()); + + typedef typename boost::property_traits::value_type Point_3; + typedef typename boost::property_traits::reference Point_ref; + + std::set degenerate_face_set; + degenerate_faces(face_range, tmesh, std::inserter(degenerate_face_set, degenerate_face_set.begin()), np); + + const std::size_t faces_size = faces(tmesh).size(); + + if(degenerate_face_set.empty()) + return true; + + if(degenerate_face_set.size() == faces_size) + { + clear(tmesh); + return true; + } + + // Sanitize the face range by adding adjacent degenerate faces + const std::size_t range_size = face_range.size(); + bool is_range_full_mesh = (range_size == faces_size); + if(!is_range_full_mesh) + { + std::list faces_to_visit(degenerate_face_set.begin(), degenerate_face_set.end()); + + while(!faces_to_visit.empty()) + { + face_descriptor fd = faces_to_visit.front(); + faces_to_visit.pop_front(); + + for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + for(halfedge_descriptor inc_hd : halfedges_around_target(hd, tmesh)) + { + face_descriptor adj_fd = face(inc_hd, tmesh); + if(adj_fd == GT::null_face() || adj_fd == fd) + continue; + + if(is_degenerate_triangle_face(adj_fd, tmesh)) + { + if(degenerate_face_set.insert(adj_fd).second) + { + // successful insertion means we did not know about this face before + faces_to_visit.push_back(adj_fd); + } + } + } + } + } + } + + // Note that there can't be any null edge incident to the degenerate faces range, + // otherwise we would have a null face incident to the face range, and that is not possible + // after the sanitization above + std::set edge_range; + for(face_descriptor fd : degenerate_face_set) + for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + edge_range.insert(edge(hd, tmesh)); + + // First remove edges of length 0 + bool all_removed = remove_degenerate_edges(edge_range, tmesh, degenerate_face_set, np); + + CGAL_assertion_code(for(face_descriptor fd : degenerate_face_set) {) + CGAL_assertion(is_degenerate_triangle_face(fd, tmesh)); + CGAL_assertion_code(}) + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + { + std::cout <<"Done with null edges.\n"; + std::ofstream output("/tmp/no_null_edges.off"); + output << std::setprecision(17) << tmesh << "\n"; + output.close(); + } +#endif + + // Then, remove triangles made of 3 collinear points + + // start by filtering out border faces + // TODO: shall we avoid doing that in case a non-manifold vertex on the boundary or if a whole component disappear? + std::set border_deg_faces; + for(face_descriptor f : degenerate_face_set) + { + halfedge_descriptor h = halfedge(f, tmesh); + for(int i=0; i<3; ++i) + { + if(is_border(opposite(h, tmesh), tmesh)) + { + border_deg_faces.insert(f); + break; + } + + h = next(h, tmesh); + } + } + + while(!border_deg_faces.empty()) + { + face_descriptor f_to_rm = *border_deg_faces.begin(); + border_deg_faces.erase(border_deg_faces.begin()); + + halfedge_descriptor h = halfedge(f_to_rm, tmesh); + for(int i=0; i<3; ++i) + { + face_descriptor f = face(opposite(h, tmesh), tmesh); + if(f!=GT::null_face()) + { + if(is_degenerate_triangle_face(f, tmesh, np)) + border_deg_faces.insert(f); + } + + h = next(h, tmesh); + } + + while(!is_border(opposite(h, tmesh), tmesh)) + { + h = next(h, tmesh); + } + + degenerate_face_set.erase(f_to_rm); + Euler::remove_face(h, tmesh); + } + + // Ignore faces with null edges + if(!all_removed) + { + std::map are_degenerate_edges; + + for(face_descriptor fd : degenerate_face_set) + { + for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + edge_descriptor ed = edge(hd, tmesh); + std::pair::iterator, bool> is_insert_successful = + are_degenerate_edges.insert(std::make_pair(ed, false)); + + bool is_degenerate = false; + if(is_insert_successful.second) + { + // did not previously exist in the map, so actually have to check if it is degenerate + if(traits.equal_3_object()(get(vpmap, target(ed, tmesh)), get(vpmap, source(ed, tmesh)))) + is_degenerate = true; + } + + is_insert_successful.first->second = is_degenerate; + + if(is_degenerate) + { + halfedge_descriptor h = halfedge(ed, tmesh); + if(!is_border(h, tmesh)) + degenerate_face_set.erase(face(h, tmesh)); + + h = opposite(h, tmesh); + if(!is_border(h, tmesh)) + degenerate_face_set.erase(face(h, tmesh)); + } + } + } + } + + // first remove degree 3 vertices that are part of a cap + // (only the vertex in the middle of the opposite edge) + // This removal does not change the shape of the mesh. + while(!degenerate_face_set.empty()) + { + std::set vertices_to_remove; + for(face_descriptor fd : degenerate_face_set) + { + for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + vertex_descriptor vd = target(hd, tmesh); + if(degree(vd, tmesh) == 3 && is_border(vd, tmesh)==GT::null_halfedge()) + { + vertices_to_remove.insert(vd); + break; + } + } + } + + for(vertex_descriptor vd : vertices_to_remove) + { + for(halfedge_descriptor hd2 : halfedges_around_target(vd, tmesh)) + degenerate_face_set.erase(face(hd2, tmesh)); + + // remove the central vertex and check if the new face is degenerated + halfedge_descriptor hd = CGAL::Euler::remove_center_vertex(halfedge(vd, tmesh), tmesh); + if(is_degenerate_triangle_face(face(hd, tmesh), tmesh, np)) + { + degenerate_face_set.insert(face(hd, tmesh)); + } + } + + if(vertices_to_remove.empty()) + break; + } + + while(!degenerate_face_set.empty()) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << "Loop on removing deg faces\n"; + + // ensure the mesh is not broken + { + std::ofstream out("/tmp/out.off"); + out << tmesh; + out.close(); + + std::vector points; + std::vector > triangles; + std::ifstream in("/tmp/out.off"); + CGAL::read_OFF(in, points, triangles); + if(!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)) + { + std::cerr << "Warning: got a polygon soup (may simply be a non-manifold vertex)!\n"; + } + } +#endif + + face_descriptor fd = *degenerate_face_set.begin(); + + // look whether an incident triangle is also degenerate + bool detect_cc_of_degenerate_triangles = false; + for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + face_descriptor adjacent_face = face(opposite(hd, tmesh), tmesh); + if(adjacent_face!=GT::null_face() && degenerate_face_set.count(adjacent_face)) + { + detect_cc_of_degenerate_triangles = true; + break; + } + } + + if(!detect_cc_of_degenerate_triangles) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " no degenerate neighbors, using a flip.\n"; +#endif + degenerate_face_set.erase(degenerate_face_set.begin()); + + // flip the longest edge of the triangle + Point_ref p1 = get(vpmap, target(halfedge(fd, tmesh), tmesh)); + Point_ref p2 = get(vpmap, target(next(halfedge(fd, tmesh), tmesh), tmesh)); + Point_ref p3 = get(vpmap, source(halfedge(fd, tmesh), tmesh)); + + CGAL_assertion(p1!=p2 && p1!=p3 && p2!=p3); + + typename Traits::Compare_distance_3 compare_distance = traits.compare_distance_3_object(); + + halfedge_descriptor edge_to_flip; + if(compare_distance(p1,p2, p1,p3) != CGAL::SMALLER) // p1p2 > p1p3 + { + if(compare_distance(p1,p2, p2,p3) != CGAL::SMALLER) // p1p2 > p2p3 + // flip p1p2 + edge_to_flip = next(halfedge(fd, tmesh), tmesh); + else + // flip p2p3 + edge_to_flip = prev(halfedge(fd, tmesh), tmesh); + } + else + { + if(compare_distance(p1,p3, p2,p3) != CGAL::SMALLER) // p1p3>p2p3 + //flip p3p1 + edge_to_flip = halfedge(fd, tmesh); + else + //flip p2p3 + edge_to_flip = prev(halfedge(fd, tmesh), tmesh); + } + + face_descriptor opposite_face=face(opposite(edge_to_flip, tmesh), tmesh); + if(opposite_face == GT::null_face()) + { + // simply remove the face + Euler::remove_face(edge_to_flip, tmesh); + } + else + { + // condition for the flip to be valid (the edge to be created do not already exists) + if(!halfedge(target(next(edge_to_flip, tmesh), tmesh), + target(next(opposite(edge_to_flip, tmesh), tmesh), tmesh), + tmesh).second) + { + Euler::flip_edge(edge_to_flip, tmesh); + } + else + { + all_removed = false; +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " WARNING: flip is not possible\n"; + // \todo Let p and q be the vertices opposite to `edge_to_flip`, and let + // r be the vertex of `edge_to_flip` that is the furthest away from + // the edge `pq`. In that case I think we should remove all the triangles + // so that the triangle pqr is in the mesh. +#endif + } + } + } + else + { + // Process a connected component of degenerate faces + // get all the faces from the connected component + // and the boundary edges + std::set cc_faces; + std::vector queue; + std::vector boundary_hedges; + std::vector inside_hedges; + queue.push_back(fd); + cc_faces.insert(fd); + + while(!queue.empty()) + { + face_descriptor top=queue.back(); + queue.pop_back(); + for(halfedge_descriptor hd : halfedges_around_face(halfedge(top, tmesh), tmesh)) + { + face_descriptor adjacent_face = face(opposite(hd, tmesh), tmesh); + if(adjacent_face==GT::null_face() || degenerate_face_set.count(adjacent_face)==0) + { + boundary_hedges.push_back(hd); + } + else + { + if(cc_faces.insert(adjacent_face).second) + queue.push_back(adjacent_face); + + if(hd < opposite(hd, tmesh)) + inside_hedges.push_back(hd); + } + } + } + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " Deal with a cc of " << cc_faces.size() << " degenerate faces.\n"; + /// dump cc_faces + { + int id = 0; + std::map vids; + for(face_descriptor f : cc_faces) + { + if(vids.insert(std::make_pair(target(halfedge(f, tmesh), tmesh), id)).second) ++id; + if(vids.insert(std::make_pair(target(next(halfedge(f, tmesh), tmesh), tmesh), id)).second) ++id; + if(vids.insert(std::make_pair(target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh), id)).second) ++id; + } + + std::ofstream output("/tmp/cc_faces.off"); + output << std::setprecision(17); + output << "OFF\n" << vids.size() << " " << cc_faces.size() << " 0\n"; + std::vector points(vids.size()); + typedef std::pair Pair_type; + for(const Pair_type& p : vids) + + points[p.second] = get(vpmap, p.first); + for(const Point_3& p : points) + output << p << "\n"; + + for(face_descriptor f : cc_faces) + { + output << "3 " + << vids[target(halfedge(f, tmesh), tmesh)] << " " + << vids[target(next(halfedge(f, tmesh), tmesh), tmesh)] << " " + << vids[target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh)] << "\n"; + } + + for(std::size_t pid=2; pid!=points.size(); ++pid) + { + CGAL_assertion(collinear(points[0], points[1], points[pid])); + } + } +#endif + + // find vertices strictly inside the cc + std::set boundary_vertices; + for(halfedge_descriptor hd : boundary_hedges) + boundary_vertices.insert(target(hd, tmesh)); + + std::set inside_vertices; + for(halfedge_descriptor hd : inside_hedges) + { + if(!boundary_vertices.count(target(hd, tmesh))) + inside_vertices.insert(target(hd, tmesh)); + if(!boundary_vertices.count(source(hd, tmesh))) + inside_vertices.insert(source(hd, tmesh)); + } + + // v-e+f = 1 for a topological disk and e = (3f+#boundary_edges)/2 + if(boundary_vertices.size()+inside_vertices.size() - + (cc_faces.size()+boundary_hedges.size())/2 != 1) + { + //cc_faces does not define a topological disk + /// \todo Find to way to handle that case +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " WARNING: Cannot remove the component of degenerate faces: not a topological disk.\n"; +#endif + + for(face_descriptor f : cc_faces) + degenerate_face_set.erase(f); + + continue; + } + + // preliminary step to check if the operation is possible + // sort the boundary points along the common supporting line + // we first need a reference point + typedef internal::Less_vertex_point Less_vertex; + std::pair< + typename std::set::iterator, + typename std::set::iterator > ref_vertices = + boost::minmax_element(boundary_vertices.begin(), + boundary_vertices.end(), + Less_vertex(traits, vpmap)); + + // and then we sort the vertices using this reference point + typedef std::set Sorted_point_set; + Sorted_point_set sorted_points; + for(vertex_descriptor v : boundary_vertices) + sorted_points.insert(get(vpmap, v)); + + CGAL_assertion(sorted_points.size()== + std::set(sorted_points.begin(), + sorted_points.end()).size()); + + CGAL_assertion(get(vpmap, *ref_vertices.first) == *sorted_points.begin()); + CGAL_assertion(get(vpmap, *ref_vertices.second) == *std::prev(sorted_points.end())); + + const Point_3& xtrm1 = *sorted_points.begin(); + const Point_3& xtrm2 = *std::prev(sorted_points.end()); + + // recover halfedges on the hole, bounded by the reference vertices + std::vector side_one, side_two; + + // look for the outgoing border halfedge of the first extreme point + for(halfedge_descriptor hd : boundary_hedges) + { + if(get(vpmap, source(hd, tmesh)) == xtrm1) + { + side_one.push_back(hd); + break; + } + } + + CGAL_assertion(side_one.size() == 1); + + bool non_monotone_border = false; + + while(get(vpmap, target(side_one.back(), tmesh)) != xtrm2) + { + vertex_descriptor prev_vertex = target(side_one.back(), tmesh); + for(halfedge_descriptor hd : boundary_hedges) + { + if(source(hd, tmesh) == prev_vertex) + { + if(get(vpmap, target(hd, tmesh)) < get(vpmap, prev_vertex)) + non_monotone_border = true; + + side_one.push_back(hd); + break; + } + } + + if(non_monotone_border) + break; + } + + if(non_monotone_border) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " WARNING: Cannot remove the component of degenerate faces: border not a monotonic cycle.\n"; +#endif + + for(face_descriptor f : cc_faces) + degenerate_face_set.erase(f); + + continue; + } + + // look for the outgoing border halfedge of second extreme vertex + for(halfedge_descriptor hd : boundary_hedges) + { + if(source(hd, tmesh) == target(side_one.back(), tmesh)) + { + side_two.push_back(hd); + break; + } + } + + CGAL_assertion(side_two.size() == 1); + + while(target(side_two.back(), tmesh) != source(side_one.front(), tmesh)) + { + vertex_descriptor prev_vertex = target(side_two.back(), tmesh); + for(halfedge_descriptor hd : boundary_hedges) + { + if(source(hd, tmesh) == prev_vertex) + { + if(get(vpmap, target(hd, tmesh)) > get(vpmap, prev_vertex)) + non_monotone_border = true; + + side_two.push_back(hd); + break; + } + } + + if(non_monotone_border) + break; + } + + if(non_monotone_border) + { +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " WARNING: Cannot remove the component of degenerate faces: border not a monotonic cycle.\n"; +#endif + + for(face_descriptor f : cc_faces) + degenerate_face_set.erase(f); + + continue; + } + + CGAL_assertion(side_one.size()+side_two.size()==boundary_hedges.size()); + + // reverse the order of the second side so as to follow + // the same order than side one + std::reverse(side_two.begin(), side_two.end()); + for(halfedge_descriptor& h : side_two) + h = opposite(h, tmesh); + + //make sure the points of the vertices along side_one are correctly sorted + std::vector side_points; + side_points.reserve(side_one.size()+1); + side_points.push_back(get(vpmap, source(side_one.front(), tmesh))); + + for(halfedge_descriptor h : side_one) + side_points.push_back(get(vpmap, target(h, tmesh))); + + CGAL_assertion(get(vpmap,source(side_one.front(), tmesh)) == side_points.front()); + CGAL_assertion(get(vpmap,target(side_one.back(), tmesh)) == side_points.back()); + + //\todo the reordering could lead to the apparition of null edges. + std::sort(side_points.begin(), side_points.end()); + + CGAL_assertion(std::unique(side_points.begin(), side_points.end())==side_points.end()); + for(std::size_t i=0; i::iterator side_one_it = side_one.begin(); + typename std::vector::iterator side_two_it = side_two.begin(); + for(;it_pt!=it_pt_end;++it_pt) + { + // check if it_pt is the point of the target of one or two halfedges + bool target_of_side_one = (get(vpmap, target(*side_one_it, tmesh)) == *it_pt); + bool target_of_side_two = (get(vpmap, target(*side_two_it, tmesh)) == *it_pt); + + if(target_of_side_one && target_of_side_two) + { + for(halfedge_descriptor h : halfedges_around_target(*side_one_it, tmesh)) + { + if(source(h, tmesh) == target(*side_two_it, tmesh)) + { + non_collapsable = true; + break; + } + } + } + else + { + CGAL_assertion(target_of_side_one || target_of_side_two); + vertex_descriptor v1 = target_of_side_one ? target(*side_one_it, tmesh) + : target(*side_two_it, tmesh); + vertex_descriptor v2 = target_of_side_two ? target(next(opposite(*side_one_it, tmesh), tmesh), tmesh) + : target(next(*side_two_it, tmesh), tmesh); + for(halfedge_descriptor h : halfedges_around_target(v1, tmesh)) + { + if(source(h, tmesh)==v2) + { + non_collapsable=true; + break; + } + } + } + + if(non_collapsable) break; + if(target_of_side_one) ++side_one_it; + if(target_of_side_two) ++side_two_it; + } + + if(non_collapsable) + { + for(face_descriptor f : cc_faces) + degenerate_face_set.erase(f); + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " WARNING: cannot remove a connected components of degenerate faces.\n"; +#endif + continue; + } + + // now proceed to the fix + // update the face and halfedge vertex pointers on the boundary + for(halfedge_descriptor h : boundary_hedges) + { + set_face(h, GT::null_face(), tmesh); + set_halfedge(target(h, tmesh), h, tmesh); + } + + // update next/prev pointers of boundary_hedges + for(halfedge_descriptor h : boundary_hedges) + { + halfedge_descriptor next_candidate = next(h, tmesh); + while(face(next_candidate, tmesh)!=GT::null_face()) + next_candidate = next(opposite(next_candidate, tmesh), tmesh); + + set_next(h, next_candidate, tmesh); + } + + // remove degenerate faces + for(face_descriptor f : cc_faces) + { + degenerate_face_set.erase(f); + remove_face(f, tmesh); + } + + // remove interior edges + for(halfedge_descriptor h : inside_hedges) + remove_edge(edge(h, tmesh), tmesh); + + // remove interior vertices + for(vertex_descriptor v : inside_vertices) + remove_vertex(v, tmesh); + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + std::cout << " side_one.size() " << side_one.size() << "\n"; + std::cout << " side_two.size() " << side_two.size() << "\n"; +#endif + + CGAL_assertion(source(side_one.front(), tmesh) == *ref_vertices.first); + CGAL_assertion(source(side_two.front(), tmesh) == *ref_vertices.first); + CGAL_assertion(target(side_one.back(), tmesh) == *ref_vertices.second); + CGAL_assertion(target(side_two.back(), tmesh) == *ref_vertices.second); + + // now split each side to contains the same sequence of points + // first side + int hi = 0; + + for(typename Sorted_point_set::iterator it=std::next(sorted_points.begin()), + it_end=sorted_points.end(); it!=it_end; ++it) + { + CGAL_assertion(*std::prev(it) == get(vpmap, source(side_one[hi], tmesh))); + + if(*it != get(vpmap, target(side_one[hi], tmesh))) + { + // split the edge and update the point + halfedge_descriptor h1 = next(opposite(side_one[hi], tmesh), tmesh); + put(vpmap, target(Euler::split_edge(side_one[hi], tmesh), tmesh), *it); + + // split_edge updates the halfedge of the source vertex of h, + // since we reuse later the halfedge of the first refernce vertex + // we must set it as we need. + if(source(h1, tmesh) == *ref_vertices.first) + set_halfedge(*ref_vertices.first, prev(prev(side_one[hi], tmesh), tmesh), tmesh); + + // retriangulate the opposite face + if(face(h1, tmesh) != GT::null_face()) + Euler::split_face(h1, opposite(side_one[hi], tmesh), tmesh); + } + else + { + ++hi; + } + } + + // second side + hi = 0; + for(typename Sorted_point_set::iterator it=std::next(sorted_points.begin()), + it_end=sorted_points.end(); it!=it_end; ++it) + { + CGAL_assertion(*std::prev(it) == get(vpmap, source(side_two[hi], tmesh))); + if(*it != get(vpmap, target(side_two[hi], tmesh))) + { + // split the edge and update the point + halfedge_descriptor h2 = Euler::split_edge(side_two[hi], tmesh); + put(vpmap, target(h2, tmesh), *it); + + // split_edge updates the halfedge of the source vertex of h, + // since we reuse later the halfedge of the first refernce vertex + // we must set it as we need. + if(source(h2, tmesh) == *ref_vertices.first) + set_halfedge(*ref_vertices.first, opposite(h2, tmesh), tmesh); + + // retriangulate the face + if(face(h2, tmesh) != GT::null_face()) + Euler::split_face(h2, next(side_two[hi], tmesh), tmesh); + } + else + { + ++hi; + } + } + + CGAL_assertion(target(halfedge(*ref_vertices.first, tmesh), tmesh) == *ref_vertices.first); + CGAL_assertion(face(halfedge(*ref_vertices.first, tmesh), tmesh) == GT::null_face()); + +#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG + { + halfedge_descriptor h_side2 = halfedge(*ref_vertices.first, tmesh); + halfedge_descriptor h_side1 = next(h_side2, tmesh); + + do + { + CGAL_assertion(get(vpmap, source(h_side1, tmesh)) == get(vpmap, target(h_side2, tmesh))); + CGAL_assertion(get(vpmap, target(h_side1, tmesh)) == get(vpmap, source(h_side2, tmesh))); + + if(target(next(opposite(h_side1, tmesh), tmesh), tmesh) == + target(next(opposite(h_side2, tmesh), tmesh), tmesh)) + { + CGAL_assertion(!"Forbidden simplification"); + } + + h_side2 = prev(h_side2, tmesh); + h_side1 = next(h_side1, tmesh); + } + while(target(h_side1, tmesh) != *ref_vertices.second); + } +#endif + + // remove side1 and replace its opposite hedges by those of side2 + halfedge_descriptor h_side2 = halfedge(*ref_vertices.first, tmesh); + halfedge_descriptor h_side1 = next(h_side2, tmesh); + for(;;) + { + CGAL_assertion(get(vpmap, source(h_side1, tmesh)) == get(vpmap, target(h_side2, tmesh))); + CGAL_assertion(get(vpmap, target(h_side1, tmesh)) == get(vpmap, source(h_side2, tmesh))); + + // backup target vertex + vertex_descriptor vertex_to_remove = target(h_side1, tmesh); + if(vertex_to_remove!=*ref_vertices.second) + { + vertex_descriptor replacement_vertex = source(h_side2, tmesh); + // replace the incident vertex + for(halfedge_descriptor hd : halfedges_around_target(h_side1, tmesh)) + set_target(hd, replacement_vertex, tmesh); + } + + // prev side2 hedge for next loop + halfedge_descriptor h_side2_for_next_turn = prev(h_side2, tmesh); + + // replace the opposite of h_side1 by h_side2 + halfedge_descriptor opposite_h_side1 = opposite(h_side1, tmesh); + face_descriptor the_face = face(opposite_h_side1, tmesh); + set_face(h_side2, the_face, tmesh); + + if(the_face!=GT::null_face()) + set_halfedge(the_face, h_side2, tmesh); + + set_next(h_side2, next(opposite_h_side1, tmesh), tmesh); + set_next(prev(opposite_h_side1, tmesh), h_side2, tmesh); + + // take the next hedges + edge_descriptor edge_to_remove = edge(h_side1, tmesh); + h_side1 = next(h_side1, tmesh); + + // now remove the extra edge + remove_edge(edge_to_remove, tmesh); + + // ... and the extra vertex if it's not the second reference + if(vertex_to_remove == *ref_vertices.second) + { + // update the halfedge pointer of the last vertex (others were already from side 2) + CGAL_assertion(target(opposite(h_side2, tmesh), tmesh) == vertex_to_remove); + set_halfedge(vertex_to_remove, opposite(h_side2, tmesh), tmesh); + break; + } + else + { + remove_vertex(vertex_to_remove , tmesh); + } + + h_side2 = h_side2_for_next_turn; + } + } + } + + return all_removed; +} + +template +bool remove_degenerate_faces(const FaceRange& face_range, + TriangleMesh& tmesh) +{ + return remove_degenerate_faces(face_range, tmesh, CGAL::parameters::all_default()); +} + +template +bool remove_degenerate_faces(TriangleMesh& tmesh, + const CGAL_PMP_NP_CLASS& np) +{ + return remove_degenerate_faces(faces(tmesh), tmesh, np); +} + +template +bool remove_degenerate_faces(TriangleMesh& tmesh) +{ + return remove_degenerate_faces(tmesh, CGAL::parameters::all_default()); +} + +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#endif // CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_self_intersections.h new file mode 100644 index 00000000000..620ca683c7a --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_self_intersections.h @@ -0,0 +1,1772 @@ +// Copyright (c) 2015-2019 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Sebastien Loriot, +// Mael Rouxel-Labbé +// +#ifndef CGAL_POLYGON_MESH_PROCESSING_REMOVE_SELF_INTERSECTIONS_H +#define CGAL_POLYGON_MESH_PROCESSING_REMOVE_SELF_INTERSECTIONS_H + +#include + +#include +#include +#include +#include // only for the preconditions +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// @tmp to debug +#define CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + +// #define CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_SMOOTHING +// #define CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_CONSTRAINTS_IN_HOLE_FILLING + +// Self-intersection removal is done by making a big-enough hole and filling it +// +// Local self-intersection removal is more subtle and only considers self-intersections +// within a connected component. It then tries to fix those by trying successively: +// - smoothing with the sharp edges in the area being constrained +// - smoothing without the sharp edges in the area being constrained +// - hole-filling with the sharp edges in the area being constrained +// - hole-filling without the sharp edges in the area being constrained +// +// The working area grows as long as we haven't been able to fix the self-intersection, +// up to a user-defined number of times. + +namespace CGAL { +namespace Polygon_mesh_processing { +namespace internal { + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + static int self_intersections_solved_by_constrained_smoothing = 0; + static int self_intersections_solved_by_unconstrained_smoothing = 0; + static int self_intersections_solved_by_constrained_hole_filling = 0; + static int self_intersections_solved_by_unconstrained_hole_filling = 0; +#endif + +template +double compute_hausdorff_distance(const TriangleMesh& old_patch, + const TriangleMesh& new_patch) +{ +#if defined(CGAL_LINKED_WITH_TBB) + typedef CGAL::Parallel_tag Tag +#else + typedef CGAL::Sequential_tag Tag; +#endif + + const double d = Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance( + new_patch, old_patch, parameters::number_of_points_on_edges(10) + .number_of_points_on_faces(100)); + return d; +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +// @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, + const bool verbose) +{ + 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(pmesh.is_valid()); // @tmp remove + CGAL_precondition(is_valid_polygon_mesh(pmesh)); // @tmp switch to '_expensive_' + + // 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 + for(const vertex_descriptor v : border_vertices) + point_to_vs[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_vs.insert(std::make_pair(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()) + { + v = add_vertex(pmesh); + } + else + { + v = vertex_stack.back(); + vertex_stack.pop_back(); + } + + put(vpm, v, p); + } + + vface.push_back(v); + } + + patch_with_vertices.push_back(vface); + } + + typedef std::pair Vertex_pair; + typedef std::map Vertex_pair_halfedge_map; + + Vertex_pair_halfedge_map halfedge_map; + + // register border halfedges + int i = 0; + for(halfedge_descriptor h : border_hedges) + { + const vertex_descriptor vs = source(h, pmesh); + const vertex_descriptor vt = target(h, pmesh); + halfedge_map.insert(std::make_pair(std::make_pair(vs, vt), h)); + + set_halfedge(target(h, pmesh), h, pmesh); // update vertex halfedge pointer + ++i; + } + + face_descriptor f = boost::graph_traits::null_face(); +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::vector new_faces; +#endif + + for(const Vertex_face& vface : patch_with_vertices) + { + if(face_stack.empty()) + { + f = add_face(pmesh); + } + else + { + f = face_stack.back(); + face_stack.pop_back(); + } + + *out++ = f; +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + new_faces.push_back(f); +#endif + + std::vector hedges; + hedges.reserve(vface.size()); + + for(std::size_t i=0, n=vface.size(); i::null_halfedge())); + halfedge_descriptor& h = it->second; + + if(success) // this halfedge is an interior halfedge + { + if(edge_stack.empty()) + { + h = halfedge(add_edge(pmesh), pmesh); + } + else + { + h = halfedge(edge_stack.back(), pmesh); + edge_stack.pop_back(); + } + + halfedge_map[std::make_pair(vj, vi)] = opposite(h, pmesh); + } + + hedges.push_back(h); + } + + CGAL_assertion(vface.size() == hedges.size()); + + // update halfedge connections + face pointers + for(int i=0, n=vface.size(); i +FaceOutputIterator replace_faces_with_patch(const std::set::face_descriptor>& face_range, + const std::vector >& patch, + PolygonMesh& pmesh, + VPM& vpm, + FaceOutputIterator out, + const bool verbose) +{ + 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; + + std::vector border_vertices; + std::set interior_vertices; + std::vector border_hedges; + std::set interior_edges; + + for(face_descriptor fh : face_range) + { + for(halfedge_descriptor h : halfedges_around_face(halfedge(fh, pmesh), pmesh)) + { + if(halfedge(target(h, pmesh), pmesh) == h) // limit the number of insertions + interior_vertices.insert(target(h, pmesh)); + } + } + + for(face_descriptor fh : face_range) + { + for(halfedge_descriptor h : halfedges_around_face(halfedge(fh, pmesh), pmesh)) + { + CGAL_assertion(!is_border(h, pmesh)); + + const edge_descriptor e = edge(h, pmesh); + const halfedge_descriptor opp_h = opposite(h, pmesh); + const face_descriptor opp_f = face(opp_h, pmesh); + + if(is_border(opp_h, pmesh) || face_range.count(opp_f) == 0) + { + vertex_descriptor v = target(h, pmesh); + interior_vertices.erase(v); + border_hedges.push_back(h); + border_vertices.push_back(v); + } + else + { + interior_edges.insert(e); + } + } + } + + return replace_faces_with_patch(border_vertices, interior_vertices, + border_hedges, interior_edges, face_range, patch, + pmesh, vpm, out, verbose); +} + +template +void replace_faces_with_patch(const std::set::face_descriptor>& faces, + const std::vector >& patch, + PolygonMesh& pmesh, + VPM& vpm, + const bool verbose = false) +{ + CGAL::Emptyset_iterator out; + replace_faces_with_patch(faces, patch, pmesh, vpm, out, verbose); +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +template +void back_up_face_range_as_point_patch(std::vector >& point_patch, + const FaceRange& face_range, + const PolygonMesh& tmesh, + const VertexPointMap& vpmap) +{ + 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(vpmap, target(h, tmesh))); + + point_patch.push_back(face_points); + } +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +template +void constrain_sharp_and_border_edges(const FaceRange& faces, + EdgeConstrainMap& eif, + const bool constrain_sharp_edges, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt) +{ + 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 GeomTraits::Vector_3 Vector; + + std::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 0 + CGAL::Polygon_mesh_processing::detect_sharp_edges_pp(faces, tmesh, 60., eif, + parameters::weak_dihedral_angle(20.)); + + // 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 + // bounded from above + const double bound = 60.; // @fixme hardcoded + const double cos_angle = std::cos(bound * CGAL_PI / 180.); + + for(const auto& ep : is_border_of_selection) + { + bool flag = ep.second; + if(constrain_sharp_edges && !flag) + { + const halfedge_descriptor h = halfedge(ep.first, tmesh); + const face_descriptor f1 = face(h, tmesh); + const face_descriptor f2 = face(opposite(h, tmesh), tmesh); + + const Vector n1 = compute_face_normal(f1, tmesh, parameters::vertex_point_map(vpmap).geom_traits(gt)); + const Vector n2 = compute_face_normal(f2, tmesh, parameters::vertex_point_map(vpmap).geom_traits(gt)); + const FT c = gt.compute_scalar_product_3_object()(n1, n2); + + // Do not mark as sharp edges with a dihedral angle that is almost `pi` because this is likely + // due to a foldness on the mesh rather than a sharp edge that we wish to preserve + // (Ideally this would be pre-treated as part of the flatness treatment) + flag = (c <= cos_angle && c >= -cos_angle); + } + + is_border_of_selection[ep.first] = flag; // @TMP ONLY NEEDED FOR CONSTRAINED EDGES OUTPUT + put(eif, ep.first, flag); + } + + // @tmp ------------------------------ + std::ofstream out("results/constrained_edges.polylines.txt"); + out << std::setprecision(17); + for(edge_descriptor e : edges(tmesh)) + if(get(eif, e)) + out << "2 " << tmesh.point(source(e, tmesh)) << " " << tmesh.point(target(e, tmesh)) << std::endl; + out.close(); + // @tmp ------------------------------ +#endif +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +template +bool remove_self_intersections_with_smoothing(std::set::face_descriptor>& face_range, + const bool constrain_sharp_edges, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + namespace CP = CGAL::parameters; + + typedef typename boost::property_traits::value_type Point; + + typedef CGAL::Face_filtered_graph Filtered_graph; + + if(verbose) + { + std::cout << " DEBUG: repair with smoothing... (constraining sharp edges: "; + std::cout << std::boolalpha << constrain_sharp_edges << ")" << std::endl; + } + + CGAL_precondition(does_self_intersect(face_range, tmesh)); + + // Keep in memory the face patch so that we can restore the mesh if smoothing is unsatisfactory + std::vector > old_patch; + back_up_face_range_as_point_patch(old_patch, face_range, tmesh, vpmap); + + // Extract the face range as a mesh because we need to compute the Hausdorff distance + const Filtered_graph ffg(tmesh, face_range); + TriangleMesh patch_mesh; + CGAL::copy_face_graph(ffg, patch_mesh, CP::vertex_point_map(vpmap)); + + // Constrain sharp and border edges + typedef CGAL::dynamic_edge_property_t Edge_property_tag; + typedef typename boost::property_map::type EIFMap; + EIFMap eif = get(Edge_property_tag(), tmesh); + + constrain_sharp_and_border_edges(face_range, eif, constrain_sharp_edges, tmesh, vpmap, gt); + + // @todo choice of number of iterations? Till convergence && max of 100? + Polygon_mesh_processing::smooth_mesh(face_range, tmesh, CP::edge_is_constrained_map(eif) + .number_of_iterations(10) + .use_safety_constraints(false)); + + // Extract the smoothed face range (again for Hausdorff computations) + const Filtered_graph ffg_post(tmesh, face_range); + TriangleMesh post_smoothing_patch_mesh; + CGAL::copy_face_graph(ffg_post, post_smoothing_patch_mesh, CP::vertex_point_map(vpmap)); + + const bool does_patch_self_intersect = does_self_intersect(face_range, tmesh, CP::vertex_point_map(vpmap)); + const double hausdorff_dist = compute_hausdorff_distance(patch_mesh, post_smoothing_patch_mesh); + + bool is_acceptable = (!does_patch_self_intersect /*&& hausdorff_dist < 1e10*/); // @fixme hardcoded + + std::cout << "does self intersect: " << does_patch_self_intersect << std::endl; + std::cout << "hausdorff dist: " << hausdorff_dist << std::endl; + std::cout << "is_acceptable: " << is_acceptable << std::endl; + + if(!is_acceptable) // restore mesh + replace_faces_with_patch(face_range, old_patch, tmesh, vpmap, verbose); +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + else + { + if(constrain_sharp_edges) + ++self_intersections_solved_by_constrained_smoothing; + else + ++self_intersections_solved_by_unconstrained_smoothing; + } +#endif + + return is_acceptable; +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +template +bool order_border_halfedge_range(std::vector::halfedge_descriptor>& hrange, + const TriangleMesh& tmesh) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + CGAL_precondition(hrange.size() > 2); + + for(std::size_t i=0; i +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) +{ + std::ofstream out(filename); + out << std::setprecision(17); + + 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)); + + 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) + { + out << pp.first << "\n"; + pp.second = unique_id++; + } + + for(const std::vector& face : point_patch) + { + out << face.size(); + for(const Point& p : face) + out << " " << unique_points_with_id.at(p); + out << "\n"; + } + + out << std::endl; + out.close(); +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +// Hole filling can be influenced by setting a third point associated to an edge on the border of the hole. +// This third point is supposed to represent how the mesh continues on the other side of the hole. +// If that edge is a border edge, there is no third point (since the opposite face is the null face). +// Similarly if the edge is an internal sharp edge, we don't really want to use the opposite face because +// there is by definition a strong discontinuity and it might thus mislead the hole filling algorithm. +// +// Rather, we construct an artifical third point that is in the same plane as the face incident to `h`, +// defined as the third point of the imaginary equilateral triangle incident to opp(h, tmesh) +template +typename boost::property_traits::value_type +construct_artificial_third_point(const typename boost::graph_traits::halfedge_descriptor h, + const TriangleMesh& tmesh, + const VertexPointMap& vpmap, + const GeomTraits& gt) +{ + typedef typename GeomTraits::FT FT; + typedef typename boost::property_traits::value_type Point; + typedef typename boost::property_traits::reference Point_ref; + typedef typename GeomTraits::Vector_3 Vector; + + const Point_ref p1 = get(vpmap, source(h, tmesh)); + const Point_ref p2 = get(vpmap, target(h, tmesh)); + const Point_ref opp_p = get(vpmap, target(next(h, tmesh), tmesh)); + + // sqrt(3)/2 to have an equilateral triangle with p1, p2, and third_point + const FT dist = 0.5 * CGAL::sqrt(3.) * CGAL::approximate_sqrt(gt.compute_squared_distance_3_object()(p1, p2)); + + const Vector ve1 = gt.construct_vector_3_object()(p1, p2); + const Vector ve2 = gt.construct_vector_3_object()(p1, opp_p); + + // gram schmidt + const FT e1e2_sp = gt.compute_scalar_product_3_object()(ve1, ve2); + Vector orthogonalized_ve2 = gt.construct_sum_of_vectors_3_object()( + ve2, gt.construct_scaled_vector_3_object()(ve1, - e1e2_sp)); + Polygon_mesh_processing::internal::normalize(orthogonalized_ve2, gt); + + const Point mid_p1p2 = gt.construct_midpoint_3_object()(p1, p2); + const Point third_p = gt.construct_translated_point_3_object()( + mid_p1p2, gt.construct_scaled_vector_3_object()(orthogonalized_ve2, -dist)); + + return third_p; +} + +template +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, + std::vector >& point_patch, + 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, + const TriangleMesh& tmesh, + VertexPointMap& /*vpmap*/, + const GeomTraits& /*gt*/, + const bool verbose) +{ + CGAL_static_assertion((std::is_same::value_type, Point>::value)); + + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + 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 patch; + if(hole_points.size() > 3) + triangulate_hole_polyline(hole_points, third_points, std::back_inserter(patch)); + else + patch.push_back(Face_indices(0, 1, 2)); // trivial hole filling + + if(patch.empty()) + { +#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 + if(verbose) + std::cout << " DEBUG: Failed to fill a hole using Delaunay search space.\n"; + + triangulate_hole_polyline(hole_points, third_points, std::back_inserter(patch), + parameters::use_delaunay_triangulation(false)); +#endif + if(patch.empty()) + { + if(verbose) + std::cout << " DEBUG: Failed to fill a hole using the whole search space.\n"; + return false; + } + } + + // @tmp ----------------------------------------------------------- + std::cout << patch.size() << " faces in the patch" << std::endl; + std::vector > to_dump; + for(const Face_indices& face : patch) + { + std::vector point_face = { hole_points[face.first], + hole_points[face.second], + hole_points[face.third] }; + to_dump.push_back(point_face); + } + + CGAL_assertion(to_dump.size() == patch.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); + // @tmp ----------------------------------------------------------- + + // make sure that the hole filling is valid, we check that no + // edge already in the mesh is present in patch. + bool non_manifold_edge_found = false; + for(const Face_indices& triangle : 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) + { + const int vi = edges[2*k], vj = edges[2*k+1]; + + // ignore boundary edges + if(vi+1 == vj || (vj == 0 && static_cast(vi) == cc_border_vertices.size()-1)) + 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) + { + non_manifold_edge_found = true; + break; + } + } + + if(non_manifold_edge_found) + break; + } + + if(non_manifold_edge_found) + { + if(verbose) + std::cout << " DEBUG: Triangulation produced is non-manifold when plugged into the mesh.\n"; + + return false; + } + + point_patch.reserve(point_patch.size() + patch.size()); + for(const Face_indices& face : patch) + { + std::vector point_face = { hole_points[face.first], + hole_points[face.second], + hole_points[face.third] }; + point_patch.push_back(point_face); + } + + if(verbose) + std::cout << " DEBUG: Found acceptable hole-filling patch.\n"; + + return true; +} + +// 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, + std::set::vertex_descriptor>& cc_interior_vertices, + std::set::edge_descriptor>& cc_interior_edges, + std::vector::value_type> >& point_patch, + const std::vector::halfedge_descriptor>& cc_border_hedges, + const std::set::face_descriptor>& cc_faces, + const TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + + typedef typename boost::property_traits::value_type Point; + + std::vector hole_points, third_points; + hole_points.reserve(cc_border_hedges.size()); + third_points.reserve(cc_border_hedges.size()); + + for(const halfedge_descriptor h : cc_border_hedges) + { + const vertex_descriptor v = source(h, tmesh); + hole_points.push_back(get(vpmap, v)); + cc_border_vertices.push_back(v); + + CGAL_assertion(!is_border(h, tmesh)); + + if(is_border_edge(h, tmesh)) + third_points.push_back(construct_artificial_third_point(h, tmesh, vpmap, gt)); + else + third_points.push_back(get(vpmap, target(next(opposite(h, tmesh), tmesh), tmesh))); + } + + CGAL_postcondition(hole_points.size() >= 3); + + return construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges, + point_patch, hole_points, third_points, cc_border_hedges, cc_faces, + tmesh, vpmap, gt, verbose); +} + +// In that 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 +bool construct_tentative_sub_hole_patch(std::vector::value_type> >& point_patch, + const std::set::face_descriptor>& sub_cc_faces, + const std::set::face_descriptor>& cc_faces, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + 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 boost::property_traits::value_type Point; + + // Collect halfedges on the boundary of the region to be selected + // (pointing inside the domain to be remeshed) + std::set internal_hedges; + std::vector cc_border_hedges; + for(const face_descriptor fd : sub_cc_faces) + { + halfedge_descriptor h = halfedge(fd, tmesh); + for(int i=0; i<3;++i) + { + if(is_border(opposite(h, tmesh), tmesh)) + { + cc_border_hedges.push_back(h); + } + else + { + const face_descriptor opp_f = face(opposite(h, tmesh), tmesh); + if(sub_cc_faces.count(opp_f) == 0) + { + cc_border_hedges.push_back(h); + if(cc_faces.count(opp_f) != 0) + internal_hedges.insert(h); + } + } + + h = next(h, tmesh); + } + } + + // Sort halfedges so that they describe the sequence of halfedges of the hole to be made + if(!order_border_halfedge_range(cc_border_hedges, tmesh)) + { + if(verbose) + std::cout << " DEBUG: More than one border in sub-hole. Not currently handled." << std::endl; + + return false; + } + + // @todo we don't care about those sets, so instead there could be a system of output iterators + // in construct_tentative_hole_patch() instead (and here would be emptyset iterators). + std::set cc_interior_vertices; + std::set cc_interior_edges; + + std::vector cc_border_vertices; + 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()); + + for(const halfedge_descriptor h : cc_border_hedges) + { + const vertex_descriptor v = source(h, tmesh); + hole_points.push_back(get(vpmap, v)); + cc_border_vertices.push_back(v); + + CGAL_assertion(!is_border(h, tmesh)); + + if(internal_hedges.count(h) == 0 && // `h` is on the border of the full CC + !is_border_edge(h, tmesh)) + { + third_points.push_back(get(vpmap, target(next(opposite(h, tmesh), tmesh), tmesh))); + } + else // `h` is on the border of the sub CC but not on the border of the full CC + { + const Point tp = construct_artificial_third_point(h, tmesh, vpmap, gt); + third_points.push_back(tp); + } + } + + return construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges, point_patch, + hole_points, third_points, cc_border_hedges, sub_cc_faces, + tmesh, vpmap, gt, verbose); +} + +// -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- + +// This function is only called when the hole is NOT subdivided into smaller holes +template +bool fill_hole(std::vector::halfedge_descriptor>& cc_border_hedges, + std::set::face_descriptor>& cc_faces, + std::set::face_descriptor>& working_face_range, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename boost::property_traits::value_type Point; + + if(verbose) + std::cout << " DEBUG: Attempting hole-filling (no constraints), " << cc_faces.size() << " faces\n"; + + if(!order_border_halfedge_range(cc_border_hedges, tmesh)) + { + CGAL_assertion(false); // we shouldn't fail to orient the boundary cycle of the complete hole + return false; + } + + std::set cc_interior_vertices; + std::set cc_interior_edges; + + std::vector cc_border_vertices; + cc_border_vertices.reserve(cc_border_hedges.size()); + + std::vector > point_patch; + + if(!construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges, point_patch, + cc_border_hedges, cc_faces, tmesh, vpmap, gt, verbose)) + { + if(verbose) + std::cout << " DEBUG: Failed to find acceptable hole patch\n"; + + return false; + } + + // Could renew the range directly within the patch replacement function + // to avoid erasing and re-adding the same face + 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, point_patch, tmesh, vpmap, + std::inserter(working_face_range, working_face_range.end()), + verbose); + + 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_postcondition(is_valid_polygon_mesh(tmesh)); + + return true; +} + +// Same function as above but border of the hole is not known +template +bool fill_hole(std::set::face_descriptor>& cc_faces, + std::set::face_descriptor>& working_face_range, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + 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(order_border_halfedge_range(cc_border_hedges, tmesh)) + return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpmap, gt, verbose); + else + return false; +} + +// Patch is not valid if: +// - we insert the same face more than once +// - insert non-manifold edges +template +bool check_point_patch_sanity(const std::vector >& point_patch) +{ + std::set > unique_faces; + std::map, int> unique_edges; + + for(const std::vector& face : point_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; + } + + return true; +} + +template +bool fill_hole_with_constraints(std::vector::halfedge_descriptor>& cc_border_hedges, + std::set::face_descriptor>& cc_faces, + std::set::face_descriptor>& working_face_range, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename boost::property_traits::value_type Point; + + typedef CGAL::Face_filtered_graph Filtered_graph; + + if(verbose) + std::cout << " DEBUG: Attempting local hole-filling with constrained sharp edges..." << std::endl; + + // If we are treating self intersections locally, first try to constrain sharp edges in the hole + typedef CGAL::dynamic_edge_property_t Edge_property_tag; + typedef typename boost::property_map::type EIFMap; + EIFMap eif = get(Edge_property_tag(), tmesh); + + constrain_sharp_and_border_edges(cc_faces, eif, true /*constrain_sharp_edges*/, tmesh, vpmap, gt); + + // Partition the hole using these constrained edges + std::set visited_faces; + std::vector > point_patch; + + int cc_counter = 0; + for(face_descriptor f : cc_faces) + { + if(!visited_faces.insert(f).second) // already visited that face + continue; + + // gather the faces of the sub-hole + std::set 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()); + + std::cout << "CC of size " << sub_cc.size() << " (total: " << cc_faces.size() << ")" << std::endl; + ++cc_counter; + + dump_cc(sub_cc, tmesh, "results/current_cc.off"); + + // The mesh is not modified, but 'point_patch' gets filled + if(!construct_tentative_sub_hole_patch(point_patch, sub_cc, cc_faces, tmesh, vpmap, gt, verbose)) + { + // 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, vpmap, gt, verbose); + } + } + + std::cout << cc_counter << " independent sub holes" << std::endl; + + // We're assembling multiple patches so we could have the same face appearing multiple times... + if(!check_point_patch_sanity(point_patch)) + { + std::cout << "unhealthy point patch" << std::endl; + return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpmap, gt, verbose); + } + + // Keep in memory the old patch in case something goes wrong + // or if the result is unsatisfactory (self-intersections, big hausdorff distance) + std::vector > old_patch; + back_up_face_range_as_point_patch(old_patch, cc_faces, tmesh, vpmap); + + // Extract the face range as a mesh because we need to compute the Hausdorff distance + const Filtered_graph ffg(tmesh, cc_faces); + TriangleMesh patch_mesh; + CGAL::copy_face_graph(ffg, patch_mesh, parameters::vertex_point_map(vpmap)); + + // Plug the hole-filling patch in the mesh + std::set new_faces; + replace_faces_with_patch(cc_faces, point_patch, tmesh, vpmap, std::inserter(new_faces, new_faces.end()), verbose); + + // Extract the plugged-in face range as a mesh (again for Hausdorff computations) + const Filtered_graph ffg_post(tmesh, new_faces); + TriangleMesh patch_mesh_post; + CGAL::copy_face_graph(ffg_post, patch_mesh_post, parameters::vertex_point_map(vpmap)); + + const bool does_patch_self_intersect = does_self_intersect(new_faces, tmesh, CGAL::parameters::vertex_point_map(vpmap)); + const double hausdorff_dist = compute_hausdorff_distance(patch_mesh, patch_mesh_post); + + const bool is_acceptable = (!does_patch_self_intersect /*&& hausdorff_dist < 1e10*/); + + for(const face_descriptor f : cc_faces) + working_face_range.erase(f); + + if(!is_acceptable) + { + if(verbose) + std::cout << " DEBUG: Failed to fill hole with constrained sharp edges. Using basic version." << std::endl; + + // 'old_faces' is equal to 'cc_faces' but needs to be grabbed again since the mesh has been modified + std::set old_faces; + replace_faces_with_patch(new_faces, old_patch, tmesh, vpmap, std::inserter(old_faces, old_faces.end()), verbose); + working_face_range.insert(old_faces.begin(), old_faces.end()); + + return fill_hole(old_faces, working_face_range, tmesh, vpmap, gt, verbose); + } + else + { + working_face_range.insert(new_faces.begin(), new_faces.end()); + return true; + } +} + +template +bool remove_self_intersections_with_hole_filling(std::vector::halfedge_descriptor>& cc_border_hedges, + std::set::face_descriptor>& cc_faces, + std::set::face_descriptor>& working_face_range, + const bool local_self_intersection_removal, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const bool verbose) +{ + // @tmp ------------------------------ + std::ofstream out("results/zone_border.polylines.txt"); + out << std::setprecision(17); + for(const auto& h : cc_border_hedges) + out << "2 " << tmesh.point(source(h, tmesh)) << " " << tmesh.point(target(h, tmesh)) << std::endl; + out.close(); + // @tmp ------------------------------ + + bool success; + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_CONSTRAINTS_IN_HOLE_FILLING + if(false) +#else + // Do not try to impose sharp edge constraints if we are not doing local-only self intersections removal + if(local_self_intersection_removal) +#endif + success = fill_hole_with_constraints(cc_border_hedges, cc_faces, working_face_range, tmesh, vpmap, gt, verbose); + else + success = fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpmap, gt, verbose); + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + if(success) + { + #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_CONSTRAINTS_IN_HOLE_FILLING + if(local_self_intersection_removal) + ++self_intersections_solved_by_constrained_hole_filling; + else + #endif + ++self_intersections_solved_by_unconstrained_hole_filling; + } +#endif + + return success; +} + +// the parameter `step` controls how many extra layers of faces we take around the range `faces_to_remove` +template +std::pair +remove_self_intersections_one_step(std::set::face_descriptor>& faces_to_remove, + std::set::face_descriptor>& working_face_range, + TriangleMesh& tmesh, + VertexPointMap& vpmap, + const GeomTraits& gt, + const int step, + const bool preserve_genus, + const bool only_treat_self_intersections_locally, + const bool verbose) +{ + 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; + + std::set faces_to_remove_copy = faces_to_remove; + + if(verbose) + { + std::cout << "DEBUG: running remove_self_intersections_one_step, step " << step + << " with " << faces_to_remove.size() << " intersecting faces\n"; + } + + CGAL_assertion(tmesh.is_valid()); + + bool something_was_done = false; // indicates if a region was successfully remeshed + bool all_fixed = true; // indicates if all removal went well + // indicates if a removal was not possible because the region handle has + // some boundary cycle of halfedges + bool topology_issue = false; + + if(verbose) + { + std::cout << " DEBUG: is_valid in one_step(tmesh)? "; + std::cout.flush(); + std::cout << is_valid_polygon_mesh(tmesh) << "\n"; + } + + while(!faces_to_remove.empty()) + { + if(verbose) + std::cout << " DEBUG: --------------- " << faces_to_remove.size() << " faces to remove (step: " << step << ")\n"; + + // Process a connected component of faces to remove. + // collect all the faces from the connected component + std::set cc_faces; + std::vector queue(1, *faces_to_remove.begin()); // temporary queue + cc_faces.insert(queue.back()); + while(!queue.empty()) + { + face_descriptor top = queue.back(); + queue.pop_back(); + halfedge_descriptor h = halfedge(top, tmesh); + for(int i=0; i<3; ++i) + { + face_descriptor adjacent_face = face(opposite(h, tmesh), tmesh); + if(adjacent_face!=boost::graph_traits::null_face()) + { + if(faces_to_remove.count(adjacent_face) != 0 && cc_faces.insert(adjacent_face).second) + queue.push_back(adjacent_face); + } + + h = next(h, tmesh); + } + } + + if(verbose) + { + std::cout << " DEBUG: " << cc_faces.size() << " faces in CC\n"; + std::cout << " DEBUG: first face: " << get(vpmap, source(halfedge(*(cc_faces.begin()), tmesh), tmesh)) << " " + << get(vpmap, target(halfedge(*(cc_faces.begin()), tmesh), tmesh)) << " " + << get(vpmap, target(next(halfedge(*(cc_faces.begin()), tmesh), tmesh), tmesh)) << "\n"; + } + + // expand the region to be filled + if(step > 0) + { + expand_face_selection(cc_faces, tmesh, step, + make_boolean_property_map(cc_faces), + Emptyset_iterator()); + } + + + // try to compactify the selection region by also selecting all the faces included + // in the bounding box of the initial selection + std::vector stack_for_expension; + Bbox_3 bb; + for(face_descriptor fd : cc_faces) + { + for(halfedge_descriptor h : halfedges_around_face(halfedge(fd, tmesh), tmesh)) + { + bb += get(vpmap, target(h, tmesh)).bbox(); + 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)); + } + } + } + + while(!stack_for_expension.empty()) + { + halfedge_descriptor h = stack_for_expension.back(); + stack_for_expension.pop_back(); + if(cc_faces.count(face(h, tmesh)) == 1) + continue; + + if(do_overlap(bb, get(vpmap, target(next(h, tmesh), tmesh)).bbox())) + { + cc_faces.insert(face(h, tmesh)); + halfedge_descriptor candidate = opposite(next(h, tmesh), tmesh); + if(face(candidate, tmesh) != boost::graph_traits::null_face()) + stack_for_expension.push_back(candidate); + + candidate = opposite(prev(h, tmesh), tmesh); + if(face(candidate, tmesh) != boost::graph_traits::null_face()) + stack_for_expension.push_back(candidate); + } + } + + if(only_treat_self_intersections_locally) + { + if(!does_self_intersect(cc_faces, tmesh, parameters::vertex_point_map(vpmap).geom_traits(gt))) + { + if(verbose) + std::cout << " DEBUG: No self-intersection in CC\n"; + + for(const face_descriptor f : cc_faces) + faces_to_remove.erase(f); + + continue; + } + } + + if(verbose) + std::cout << " DEBUG: " << cc_faces.size() << " faces in expanded CC\n"; + + // remove faces from the set to process + for(const face_descriptor f : cc_faces) + faces_to_remove.erase(f); + + // @tmp ----------------------------------------------------- + std::stringstream ex_oss; + std::cout << "Output FULLY expanded CC #" << exp_cc_id << std::endl; + ex_oss << "results/fully_expanded_cc_" << exp_cc_id++ << ".off" << std::ends; + dump_cc(cc_faces, tmesh, ex_oss.str().c_str()); + // @tmp ----------------------------------------------------- + + //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); + + h = next(h, tmesh); + } + } + + // 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; + if(verbose) + std::cout << " DEBUG: CC not handled due to the presence at least one non-manifold vertex\n"; + + 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 + continue; + + working_face_range.insert(cc_faces.begin(), cc_faces.end()); + + // Now, we have a proper selection that we can work on. + +#ifndef CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_SMOOTHING + // First, try to smooth if we only care about local self-intersections + // Two different approaches: + // - First, try to constrain edges that are in the zone to smooth and whose dihedral angle is large, + // but not too large (we don't want to constrain edges that are created by foldings) + // - If that fails, try to smooth without any constraints, but make sure that the deviation from + // the first zone is small + // + // If smoothing fails, the face patch is restored to its pre-smoothing state. + // + // Note that there is no need to update the working range because smoothing doesn`t change + // the number of faces (and old faces are re-used). + bool fixed_by_smoothing = false; + + if(only_treat_self_intersections_locally) + { + fixed_by_smoothing = remove_self_intersections_with_smoothing(cc_faces, true, tmesh, vpmap, gt, verbose); + if(!fixed_by_smoothing) // try again, but without constraining sharp edges + { + if(verbose) + std::cout << " DEBUG: Could not be solved via smoothing with constraints\n"; + + fixed_by_smoothing = remove_self_intersections_with_smoothing(cc_faces, false, tmesh, vpmap, gt, verbose); + } + } + + if(fixed_by_smoothing) + { + if(verbose) + std::cout << " DEBUG: Solved with smoothing!\n"; + + something_was_done = true; + continue; + } + else if(verbose) + { + std::cout << " DEBUG: Could not be solved via smoothing\n"; + } +#endif + + if(verbose) + std::cout << " DEBUG: Trying hole-filling based approach\n"; + + if(!is_selection_a_topological_disk(cc_faces, tmesh)) + { + // check if the selection contains cycles of border halfedges + bool only_border_edges = true; + std::set mesh_border_hedge; + + for(halfedge_descriptor h : cc_border_hedges) + { + if(!is_border(opposite(h, tmesh), tmesh)) + only_border_edges = false; + else + mesh_border_hedge.insert(opposite(h, tmesh)); + } + + int nb_cycles = 0; + 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 + ++nb_cycles; + break; + } + else + { + typename std::set::iterator it = mesh_border_hedge.find(h); + if(it == mesh_border_hedge.end()) + break; // not a cycle + + mesh_border_hedge.erase(it); + } + } + while(true); + } + + if(nb_cycles > (only_border_edges ? 1 : 0)) + { + if(verbose) + std::cout << " DEBUG: CC not handled due to the presence of " + << nb_cycles << " of boundary edges\n"; + + topology_issue = true; + continue; + } + else + { + if(preserve_genus) + { + if(verbose) + std::cout << " DEBUG: CC not handled because it is not a topological disk (preserve_genus=true)\n"; + + all_fixed = false; + continue; + } + + // 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); + + 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) + { + if(verbose) + std::cout << " DEBUG: CC not handled because it is not a topological disk(" + << nbc << " boundary cycles)\n"; + + all_fixed = false; + continue; + } + else + { + if(verbose) + std::cout << " DEBUG: CC that is not a topological disk but has only one boundary cycle(preserve_genus=false)\n"; + } + } + } + + // sort halfedges so that they describe the sequence of halfedges of the hole to be made + if(!remove_self_intersections_with_hole_filling(cc_border_hedges, cc_faces, working_face_range, + only_treat_self_intersections_locally, + tmesh, vpmap, gt, verbose)) + { + if(verbose) + std::cout << " DEBUG: Failed to fill hole\n"; + + all_fixed = false; + } + else + { + something_was_done = true; + } + } + + if(!something_was_done) + { + faces_to_remove.swap(faces_to_remove_copy); + if(verbose) + std::cout << " DEBUG: Nothing was changed during this step, self-intersections won`t be recomputed." << std::endl; + } + + std::stringstream oss; + oss << "results/after_step_" << step << ".off" << std::ends; + std::ofstream(oss.str().c_str()) << std::setprecision(17) << tmesh; + + return std::make_pair(all_fixed, topology_issue); +} + +} // namespace internal + +/// \cond SKIP_IN_MANUAL + +template +bool remove_self_intersections(const FaceRange& face_range, + TriangleMesh& tmesh, + const NamedParameters& np) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef boost::graph_traits graph_traits; + typedef typename graph_traits::face_descriptor face_descriptor; + + // named parameter extraction + typedef typename GetVertexPointMap::type VertexPointMap; + VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, tmesh)); + + typedef typename GetGeomTraits::type GeomTraits; + GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GeomTraits()); + + bool verbose = choose_parameter(get_parameter(np, internal_np::verbosity_level), 0) > 0; + bool preserve_genus = choose_parameter(get_parameter(np, internal_np::preserve_genus), true); + + const bool only_treat_self_intersections_locally = true; + verbose = true; // @tmp + + // When treating intersections locally, we don't want to grow the working range too much + const int default_max_step = only_treat_self_intersections_locally ? 2 : 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'? + const double strong_dihedral_angle = choose_parameter(get_parameter(np, internal_np::with_dihedral_angle), 60.); + const double weak_dihedral_angle = choose_parameter(get_parameter(np, internal_np::weak_dihedral_angle), 20.); + + std::set working_face_range(face_range.begin(), face_range.end()); + + if(verbose) + std::cout << "DEBUG: Starting remove_self_intersections, is_valid(tmesh)? " << is_valid_polygon_mesh(tmesh) << "\n"; + + if(!preserve_genus) + duplicate_non_manifold_vertices(tmesh, np); + + // Look for self-intersections in the mesh and remove them + 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; + + while(++step < max_steps) + { + if(faces_to_remove.empty()) // the previous round might have been blocked due to topological constraints + { + typedef std::pair Face_pair; + std::vector self_inter; + // TODO : possible optimization to reduce the range to check with the bbox + // of the previous patches or something. + self_intersections(working_face_range, tmesh, std::back_inserter(self_inter)); + + std::cout << self_inter.size() << " intersecting pairs" << std::endl; + + for(const Face_pair& fp : self_inter) + { + faces_to_remove.insert(fp.first); + faces_to_remove.insert(fp.second); + } + } + + if(faces_to_remove.empty() && all_fixed) + { + if(verbose) + std::cout << "DEBUG: There are no more faces to remove." << std::endl; + break; + } + + std::tie(all_fixed, topology_issue) = + internal::remove_self_intersections_one_step( + faces_to_remove, working_face_range, tmesh, vpm, gt, + step, preserve_genus, only_treat_self_intersections_locally, verbose); + + if(all_fixed && topology_issue) + { + if(verbose) + std::cout << "DEBUG: boundary cycles of boundary edges involved in self-intersections.\n"; + } + } + +#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG + std::cout << "solved by constrained smoothing: " << internal::self_intersections_solved_by_constrained_smoothing << std::endl; + 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; +#endif + + return step < max_steps; +} + +template +bool remove_self_intersections(const FaceRange& face_range, TriangleMesh& tmesh) +{ + return remove_self_intersections(face_range, tmesh, parameters::all_default()); +} + +template +bool remove_self_intersections(TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np) +{ + return remove_self_intersections(faces(tmesh), tmesh, np); +} + +template +bool remove_self_intersections(TriangleMesh& tmesh) +{ + return remove_self_intersections(tmesh, parameters::all_default()); +} + +/// \endcond + +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#endif // CGAL_POLYGON_MESH_PROCESSING_REMOVE_SELF_INTERSECTIONS_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h index f635af7bc95..c73c44048c2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h @@ -16,50 +16,16 @@ #include -#include -#include -#include -#include -#include -#include - -// headers for self-intersection removal -#include -#include -#include - -#include +#include +#include +#include #include -#include -#include #include -#include -#include +#include -#include - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG -#include -#include -#endif - -#include -#include -#include - -#include -#include -#include -#include -#include #include -#ifdef DOXYGEN_RUNNING -#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters -#define CGAL_PMP_NP_CLASS NamedParameters -#endif - namespace CGAL { namespace Polygon_mesh_processing { @@ -340,2899 +306,7 @@ std::size_t remove_connected_components_of_negligible_size(TriangleMesh& tmesh) return remove_connected_components_of_negligible_size(tmesh, parameters::all_default()); } -namespace debug{ - - template - std::ostream& dump_edge_neighborhood( - typename boost::graph_traits::edge_descriptor ed, - TriangleMesh& tmesh, - const VertexPointMap& vpmap, - std::ostream& out) - { - typedef boost::graph_traits GT; - typedef typename GT::halfedge_descriptor halfedge_descriptor; - typedef typename GT::vertex_descriptor vertex_descriptor; - typedef typename GT::face_descriptor face_descriptor; - - halfedge_descriptor h = halfedge(ed, tmesh); - - std::map vertices; - std::set faces; - int vindex=0; - for(halfedge_descriptor hd : halfedges_around_target(h, tmesh)) - { - if ( vertices.insert(std::make_pair(source(hd, tmesh), vindex)).second ) - ++vindex; - if (!is_border(hd, tmesh)) - faces.insert( face(hd, tmesh) ); - } - - h=opposite(h, tmesh); - for(halfedge_descriptor hd : halfedges_around_target(h, tmesh)) - { - if ( vertices.insert(std::make_pair(source(hd, tmesh), vindex)).second ) - ++vindex; - if (!is_border(hd, tmesh)) - faces.insert( face(hd, tmesh) ); - } - - std::vector ordered_vertices(vertices.size()); - typedef std::pair Pair_type; - for(const Pair_type& p : vertices) - ordered_vertices[p.second]=p.first; - - out << "OFF\n" << ordered_vertices.size() << " " << faces.size() << " 0\n"; - for(vertex_descriptor vd : ordered_vertices) - out << get(vpmap, vd) << "\n"; - for(face_descriptor fd : faces) - { - out << "3"; - h=halfedge(fd,tmesh); - for(halfedge_descriptor hd : halfedges_around_face(h, tmesh)) - out << " " << vertices[target(hd, tmesh)]; - out << "\n"; - } - return out; - } - - template - void dump_cc_faces(const FaceRange& cc_faces, const TriangleMesh& tm, std::ostream& output) - { - typedef typename boost::property_map::const_type Vpm; - typedef typename boost::property_traits::value_type Point_3; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - Vpm vpm = get(boost::vertex_point, tm); - - int id=0; - std::map vids; - for(face_descriptor f : cc_faces) - { - if ( vids.insert( std::make_pair( target(halfedge(f, tm), tm), id) ).second ) ++id; - if ( vids.insert( std::make_pair( target(next(halfedge(f, tm), tm), tm), id) ).second ) ++id; - if ( vids.insert( std::make_pair( target(next(next(halfedge(f, tm), tm), tm), tm), id) ).second ) ++id; - } - output << std::setprecision(17); - output << "OFF\n" << vids.size() << " " << cc_faces.size() << " 0\n"; - std::vector points(vids.size()); - typedef std::pair Pair_type; - for(Pair_type p : vids) - points[p.second]=get(vpm, p.first); - for(Point_3 p : points) - output << p << "\n"; - for(face_descriptor f : cc_faces) - { - output << "3 " - << vids[ target(halfedge(f, tm), tm) ] << " " - << vids[ target(next(halfedge(f, tm), tm), tm) ] << " " - << vids[ target(next(next(halfedge(f, tm), tm), tm), tm) ] << "\n"; - } - } - -} //end of namespace debug - -namespace internal { - -template -struct Less_vertex_point{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - const Traits& m_traits; - const VertexPointMap& m_vpmap; - Less_vertex_point(const Traits& traits, const VertexPointMap& vpmap) - : m_traits(traits) - , m_vpmap(vpmap) {} - bool operator()(vertex_descriptor v1, vertex_descriptor v2) const{ - return m_traits.less_xyz_3_object()(get(m_vpmap, v1), get(m_vpmap, v2)); - } -}; - -} // end namespace internal - -/// \ingroup PMP_repairing_grp -/// collects the degenerate edges within a given range of edges. -/// -/// @tparam EdgeRange a model of `Range` with value type `boost::graph_traits::%edge_descriptor` -/// @tparam TriangleMesh a model of `HalfedgeGraph` -/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" -/// -/// @param edges a subset of edges of `tm` -/// @param tm a triangle mesh -/// @param out an output iterator in which the degenerate edges are written -/// @param np optional \ref pmp_namedparameters "Named Parameters" described below -/// -/// \cgalNamedParamsBegin -/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`. -/// The type of this map is model of `ReadWritePropertyMap`. -/// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh` -/// \cgalParamEnd -/// \cgalParamBegin{geom_traits} a geometric traits class instance. -/// The traits class must provide the nested type `Point_3`, -/// and the nested functor `Equal_3` to check whether two points are identical. -/// \cgalParamEnd -/// \cgalNamedParamsEnd -template -OutputIterator degenerate_edges(const EdgeRange& edges, - const TriangleMesh& tm, - OutputIterator out, - const NamedParameters& np) -{ - typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - - for(edge_descriptor ed : edges) - { - if(is_degenerate_edge(ed, tm, np)) - *out++ = ed; - } - return out; -} - -template -OutputIterator degenerate_edges(const EdgeRange& edges, - const TriangleMesh& tm, - OutputIterator out, - typename boost::enable_if< - typename boost::has_range_iterator - >::type* = 0) -{ - return degenerate_edges(edges, tm, out, CGAL::parameters::all_default()); -} - -/// \ingroup PMP_repairing_grp -/// calls the function `degenerate_edges()` with the range: `edges(tm)`. -/// -/// See above for the comprehensive description of the parameters. -/// -template -OutputIterator degenerate_edges(const TriangleMesh& tm, - OutputIterator out, - const NamedParameters& np -#ifndef DOXYGEN_RUNNING - , typename boost::disable_if< - boost::has_range_iterator - >::type* = 0 -#endif - ) -{ - return degenerate_edges(edges(tm), tm, out, np); -} - -template -OutputIterator -degenerate_edges(const TriangleMesh& tm, OutputIterator out) -{ - return degenerate_edges(edges(tm), tm, out, CGAL::parameters::all_default()); -} - -/// \ingroup PMP_repairing_grp -/// collects the degenerate faces within a given range of faces. -/// -/// @tparam FaceRange a model of `Range` with value type `boost::graph_traits::%face_descriptor` -/// @tparam TriangleMesh a model of `FaceGraph` -/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" -/// -/// @param faces a subset of faces of `tm` -/// @param tm a triangle mesh -/// @param out an output iterator in which the degenerate faces are put -/// @param np optional \ref pmp_namedparameters "Named Parameters" described below -/// -/// \cgalNamedParamsBegin -/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`. -/// The type of this map is model of `ReadWritePropertyMap`. -/// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh` -/// \cgalParamEnd -/// \cgalParamBegin{geom_traits} a geometric traits class instance. -/// The traits class must provide the nested functor `Collinear_3` -/// to check whether three points are collinear. -/// \cgalParamEnd -/// \cgalNamedParamsEnd -/// -template -OutputIterator degenerate_faces(const FaceRange& faces, - const TriangleMesh& tm, - OutputIterator out, - const NamedParameters& np) -{ - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - for(face_descriptor fd : faces) - { - if(is_degenerate_triangle_face(fd, tm, np)) - *out++ = fd; - } - return out; -} - -template -OutputIterator degenerate_faces(const FaceRange& faces, - const TriangleMesh& tm, - OutputIterator out, - typename boost::enable_if< - boost::has_range_iterator - >::type* = 0) -{ - return degenerate_faces(faces, tm, out, CGAL::parameters::all_default()); -} - -/// \ingroup PMP_repairing_grp -/// calls the function `degenerate_faces()` with the range: `faces(tm)`. -/// -/// See above for the comprehensive description of the parameters. -/// -template -OutputIterator degenerate_faces(const TriangleMesh& tm, - OutputIterator out, - const NamedParameters& np -#ifndef DOXYGEN_RUNNING - , typename boost::disable_if< - boost::has_range_iterator - >::type* = 0 -#endif - ) -{ - return degenerate_faces(faces(tm), tm, out, np); -} - -template -OutputIterator degenerate_faces(const TriangleMesh& tm, OutputIterator out) -{ - return degenerate_faces(faces(tm), tm, out, CGAL::parameters::all_default()); -} - -// this function removes a border edge even if it does not satisfy the link condition. -// null_vertex() is returned if the removal changes the topology of the input -template -typename boost::graph_traits::vertex_descriptor -remove_a_border_edge(typename boost::graph_traits::edge_descriptor ed, - TriangleMesh& tm, - EdgeSet& edge_set, - FaceSet& face_set) -{ - typedef boost::graph_traits GT; - typedef typename GT::edge_descriptor edge_descriptor; - typedef typename GT::halfedge_descriptor halfedge_descriptor; - typedef typename GT::face_descriptor face_descriptor; - typedef typename GT::vertex_descriptor vertex_descriptor; - - halfedge_descriptor h=halfedge(ed,tm); - - if ( is_border(h,tm) ) - h=opposite(h,tm); - - halfedge_descriptor opp_h = opposite(h,tm); - CGAL_assertion(is_border(opp_h,tm)); - CGAL_assertion(!is_border(h,tm)); - - CGAL_assertion(next(next(opp_h, tm), tm) !=opp_h); // not working for a hole made of 2 edges - CGAL_assertion(next(next(next(opp_h, tm), tm), tm) !=opp_h); // not working for a hole make of 3 edges - - if (CGAL::Euler::does_satisfy_link_condition(edge(h,tm),tm)) - { - edge_set.erase(ed); - halfedge_descriptor h=halfedge(ed, tm); - if ( is_border(h, tm) ) - h = opposite(h, tm); - - edge_set.erase(edge(prev(h, tm), tm)); - face_set.erase(face(h, tm)); - - return CGAL::Euler::collapse_edge(ed, tm); - } - - // collect edges that have one vertex in the link of - // the vertices of h and one of the vertex of h as other vertex - std::set common_incident_edges; - for(halfedge_descriptor hos : halfedges_around_source(h, tm)) - for(halfedge_descriptor hot : halfedges_around_target(h, tm)) - { - if( target(hos, tm) == source(hot, tm) ) - { - common_incident_edges.insert( edge(hot, tm) ); - common_incident_edges.insert( edge(hos, tm) ); - } - } - - CGAL_assertion(common_incident_edges.size() >=2 ); - - // in the following loop, we visit define a connected component of - // faces bounded by edges in common_incident_edges and h. We look - // for the maximal one. This set of faces is the one that will - // disappear while collapsing ed - std::set marked_faces; - - std::vector queue; - queue.push_back( opposite(next(h,tm), tm) ); - queue.push_back( opposite(prev(h,tm), tm) ); - marked_faces.insert( face(h, tm) ); - - do - { - std::vector boundary; - while(!queue.empty()) - { - halfedge_descriptor back=queue.back(); - queue.pop_back(); - face_descriptor fback=face(back,tm); - if (common_incident_edges.count(edge(back,tm))) - { - boundary.push_back(back); - continue; - } - - if ( fback==GT::null_face() || !marked_faces.insert(fback).second ) - continue; - - queue.push_back( opposite(next(back,tm), tm) ); - if ( is_border(queue.back(), tm) ) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Boundary reached during exploration, the region to be removed is not a topological disk, not handled for now.\n"; -#endif - return GT::null_vertex(); - } - - queue.push_back( opposite(prev(back,tm), tm) ); - if ( is_border(queue.back(), tm) ) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Boundary reached during exploration, the region to be removed is not a topological disk, not handled for now.\n"; -#endif - return GT::null_vertex(); - } - } - - CGAL_assertion( boundary.size() == 2 ); - common_incident_edges.erase( edge(boundary[0], tm) ); - common_incident_edges.erase( edge(boundary[1], tm) ); - if (!is_border(boundary[0], tm) || common_incident_edges.empty()) - queue.push_back(boundary[0]); - if (!is_border(boundary[1], tm) || common_incident_edges.empty()) - queue.push_back(boundary[1]); - } - while(!common_incident_edges.empty()); - - // hk1 and hk2 are bounding the region that will be removed. - // The edge of hk2 will be removed and hk2 will be replaced - // by the opposite edge of hk1 - halfedge_descriptor hk1=queue.front(); - halfedge_descriptor hk2=queue.back(); - if ( target(hk1,tm)!=source(hk2,tm) ) - std::swap(hk1, hk2); - - CGAL_assertion( target(hk1,tm)==source(hk2,tm) ); - CGAL_assertion( source(hk1,tm)==source(h,tm) ); - CGAL_assertion( target(hk2,tm)==target(h,tm) ); - - CGAL_assertion(is_valid_polygon_mesh(tm)); - if (!is_selection_a_topological_disk(marked_faces, tm)) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "The region to be removed is not a topological disk, not handled for now.\n"; -#endif - return GT::null_vertex(); - - } - - if (is_border(hk1, tm) && is_border(hk2, tm)) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "The region to be removed is an isolated region, not handled for now.\n"; -#endif - return GT::null_vertex(); - } - - // collect vertices and edges to remove and do remove faces - std::set edges_to_remove; - std::set vertices_to_remove; - for(face_descriptor fd : marked_faces) - { - halfedge_descriptor hd=halfedge(fd, tm); - for(int i=0; i<3; ++i) - { - edges_to_remove.insert( edge(hd, tm) ); - vertices_to_remove.insert( target(hd,tm) ); - hd=next(hd, tm); - } - } - - vertex_descriptor vkept=source(hk1,tm); - - //back-up next, prev halfedge pointers to be restored after removal - halfedge_descriptor hp=prev(opp_h, tm); - halfedge_descriptor hn=next(opp_h, tm); - halfedge_descriptor hk1_opp_next = next(hk2, tm); - halfedge_descriptor hk1_opp_prev = prev(hk2, tm); - face_descriptor hk1_opp_face = face(hk2,tm); - - // we will remove the target of hk2, update vertex pointers - for(halfedge_descriptor hot : halfedges_around_target(hk2, tm)) - set_target(hot, vkept, tm); - - // update halfedge pointers since hk2 will be removed - set_halfedge(vkept, opposite(hk1, tm), tm); - set_halfedge(target(hk1,tm), hk1, tm); - - // do not remove hk1 and its vertices - vertices_to_remove.erase( vkept ); - vertices_to_remove.erase( target(hk1, tm) ); - edges_to_remove.erase( edge(hk1,tm) ); - - bool hk2_equals_hp = hk2==hp; - CGAL_assertion( is_border(hk2, tm) == hk2_equals_hp ); - - /* - - case hk2!=hp - - /\ / - hk1/ \hk2 / - / \ / - ____/______\/____ - hn h_opp hp - - - case hk2==hp - - /\ - hk1/ \hk2 == hp - / \ - ____/______\ - hn h_opp - */ - - // remove vertices - for(vertex_descriptor vd : vertices_to_remove) - remove_vertex(vd, tm); - - // remove edges - for(edge_descriptor ed : edges_to_remove) - { - edge_set.erase(ed); - remove_edge(ed, tm); - } - - // remove faces - for(face_descriptor fd : marked_faces) - { - face_set.erase(fd); - remove_face(fd, tm); - } - - // now update pointers - set_face(opposite(hk1, tm), hk1_opp_face, tm); - if (!hk2_equals_hp) - { - set_next(hp, hn, tm); - set_next(opposite(hk1, tm), hk1_opp_next, tm); - set_next(hk1_opp_prev, opposite(hk1, tm), tm); - set_halfedge(hk1_opp_face, opposite(hk1, tm), tm); - } - else - { - set_next(hk1_opp_prev, opposite(hk1, tm), tm); - set_next(opposite(hk1, tm), hn, tm); - } - - CGAL_assertion(is_valid_polygon_mesh(tm)); - return vkept; -} - -template -typename boost::graph_traits::vertex_descriptor -remove_a_border_edge(typename boost::graph_traits::edge_descriptor ed, - TriangleMesh& tm) -{ - std::set::edge_descriptor> edge_set; - std::set::face_descriptor> face_set; - - return remove_a_border_edge(ed, tm, edge_set, face_set); -} - -template -bool remove_degenerate_edges(const EdgeRange& edge_range, - TriangleMesh& tmesh, - FaceSet& face_set, - const NamedParameters& np) -{ - CGAL_assertion(CGAL::is_triangle_mesh(tmesh)); - CGAL_assertion(CGAL::is_valid_polygon_mesh(tmesh)); - - using parameters::get_parameter; - using parameters::choose_parameter; - - typedef TriangleMesh TM; - typedef typename boost::graph_traits GT; - typedef typename GT::edge_descriptor edge_descriptor; - typedef typename GT::halfedge_descriptor halfedge_descriptor; - typedef typename GT::face_descriptor face_descriptor; - typedef typename GT::vertex_descriptor vertex_descriptor; - - typedef typename GetVertexPointMap::type VertexPointMap; - VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_property_map(vertex_point, tmesh)); - - typedef typename GetGeomTraits::type Traits; - - std::size_t nb_deg_faces = 0; - bool all_removed=false; - bool some_removed=true; - - bool preserve_genus = parameters::choose_parameter(parameters::get_parameter(np, internal_np::preserve_genus), true); - - // collect edges of length 0 - while(some_removed && !all_removed) - { - some_removed=false; - all_removed=true; - std::set degenerate_edges_to_remove; - degenerate_edges(edge_range, tmesh, std::inserter(degenerate_edges_to_remove, - degenerate_edges_to_remove.end())); - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Found " << degenerate_edges_to_remove.size() << " null edges.\n"; -#endif - - // first try to remove all collapsable edges - typename std::set::iterator it = degenerate_edges_to_remove.begin(); - while (it!=degenerate_edges_to_remove.end()) - { - edge_descriptor e = *it; - if (CGAL::Euler::does_satisfy_link_condition(e,tmesh)) - { - halfedge_descriptor h = halfedge(e, tmesh); - degenerate_edges_to_remove.erase(it); - - // remove edges that could also be set for removal - if ( face(h, tmesh)!=GT::null_face() ) - { - ++nb_deg_faces; - degenerate_edges_to_remove.erase(edge(prev(h, tmesh), tmesh)); - face_set.erase(face(h, tmesh)); - } - - if (face(opposite(h, tmesh), tmesh)!=GT::null_face()) - { - ++nb_deg_faces; - degenerate_edges_to_remove.erase(edge(prev(opposite(h, tmesh), tmesh), tmesh)); - face_set.erase(face(opposite(h, tmesh), tmesh)); - } - - //now remove the edge - CGAL::Euler::collapse_edge(e, tmesh); - - // some_removed is not updated on purpose because if nothing - // happens below then nothing can be done - it = degenerate_edges_to_remove.begin(); - } - else - ++it; - } - - CGAL_assertion( is_valid_polygon_mesh(tmesh) ); -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Remaining " << degenerate_edges_to_remove.size() << " null edges to be handled.\n"; -#endif - - while (!degenerate_edges_to_remove.empty()) - { - edge_descriptor ed = *degenerate_edges_to_remove.begin(); - degenerate_edges_to_remove.erase(degenerate_edges_to_remove.begin()); - - halfedge_descriptor h = halfedge(ed, tmesh); - - if (CGAL::Euler::does_satisfy_link_condition(ed,tmesh)) - { - // remove edges that could also be set for removal - if ( face(h, tmesh)!=GT::null_face() ) - { - ++nb_deg_faces; - degenerate_edges_to_remove.erase(edge(prev(h, tmesh), tmesh)); - face_set.erase(face(h, tmesh)); - } - - if (face(opposite(h, tmesh), tmesh)!=GT::null_face()) - { - ++nb_deg_faces; - degenerate_edges_to_remove.erase(edge(prev(opposite(h, tmesh), tmesh), tmesh)); - face_set.erase(face(opposite(h, tmesh), tmesh)); - } - - //now remove the edge - CGAL::Euler::collapse_edge(ed, tmesh); - some_removed = true; - } - else - { - //handle the case when the edge is incident to a triangle hole - //we first fill the hole and try again - if ( is_border(ed, tmesh) ) - { - halfedge_descriptor hd = halfedge(ed,tmesh); - if (!is_border(hd,tmesh)) - hd=opposite(hd,tmesh); - - if (is_triangle(hd, tmesh)) - { - if (!preserve_genus) - { - Euler::fill_hole(hd, tmesh); - degenerate_edges_to_remove.insert(ed); // reinsert the edge for future processing - } - else - { - all_removed=false; - } - continue; - } - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Calling remove_a_border_edge\n"; -#endif - - vertex_descriptor vd = remove_a_border_edge(ed, tmesh, degenerate_edges_to_remove, face_set); - if (vd == GT::null_vertex()) - { - // TODO: if some border edges are later removed, the edge might be processable later - // for example if it belongs to boundary cycle of edges where the number of non-degenerate - // edges is 2. That's what happen with fused_vertices.off in the testsuite where the edges - // are not processed the same way with Polyhedron and Surface_mesh. In the case of Polyhedron - // more degenerate edges could be removed. - all_removed=false; - } - else - some_removed=true; - - continue; - } - else - { - halfedge_descriptor hd = halfedge(ed,tmesh); - // if both vertices are boundary vertices we can't do anything - bool impossible = false; - for(halfedge_descriptor h : halfedges_around_target(hd, tmesh)) - { - if (is_border(h, tmesh)) - { - impossible = true; - break; - } - } - - if (impossible) - { - impossible=false; - for(halfedge_descriptor h : halfedges_around_source(hd, tmesh)) - { - if (is_border(h, tmesh)) - { - impossible = true; - break; - } - } - - if (impossible) - { - all_removed=false; - continue; - } - } - } - - // When the edge does not satisfy the link condition, it means that it cannot be - // collapsed as is. In the following if there is a topological issue - // with contracting the edge (component or geometric feature that disappears), - // nothing is done. - // We start by marking the faces that are incident to an edge endpoint. - // If the set of marked faces is a topologically disk, then we simply remove all the simplicies - // inside the disk and star the hole with the edge vertex kept. - // If the set of marked faces is not a topological disk, it has some non-manifold vertices - // on its boundary. We need to mark additional faces to make it a topological disk. - // We can then apply the star hole procedure. - // Right now we additionally mark the smallest connected components of non-marked faces - // (using the number of faces) - - //backup central point - typename Traits::Point_3 pt = get(vpmap, source(ed, tmesh)); - - // mark faces of the link of each endpoints of the edge which collapse is not topologically valid - std::set marked_faces; - - // first endpoint - for(halfedge_descriptor hd : CGAL::halfedges_around_target(halfedge(ed,tmesh), tmesh) ) - if (!is_border(hd,tmesh)) marked_faces.insert( face(hd, tmesh) ); - - // second endpoint - for(halfedge_descriptor hd : CGAL::halfedges_around_target(opposite(halfedge(ed, tmesh), tmesh), tmesh) ) - if (!is_border(hd,tmesh)) marked_faces.insert( face(hd, tmesh) ); - - // extract the halfedges on the boundary of the marked region - std::vector border; - for(face_descriptor fd : marked_faces) - { - for(halfedge_descriptor hd : CGAL::halfedges_around_face(halfedge(fd,tmesh), tmesh)) - { - halfedge_descriptor hd_opp = opposite(hd, tmesh); - if ( is_border(hd_opp, tmesh) || - marked_faces.count( face(hd_opp, tmesh) ) == 0 ) - { - border.push_back( hd ); - } - } - } - - if (border.empty() ) - { - // a whole connected component (without boundary) got selected and will disappear (not handled for now) -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Trying to remove a whole connected component, not handled yet\n"; -#endif - all_removed=false; - continue; - } - - // define cc of border halfedges: two halfedges are in the same cc - // if they are on the border of the cc of non-marked faces. - typedef CGAL::Union_find UF_ds; - UF_ds uf; - std::map handles; - // one cc per border halfedge - for(halfedge_descriptor hd : border) - handles.insert( std::make_pair(hd, uf.make_set(hd)) ); - - // join cc's - for(halfedge_descriptor hd : border) - { - CGAL_assertion( marked_faces.count( face( hd, tmesh) ) > 0); - CGAL_assertion( marked_faces.count( face( opposite(hd, tmesh), tmesh) ) == 0 ); - halfedge_descriptor candidate = hd; - - do { - candidate = prev( opposite(candidate, tmesh), tmesh ); - } while( !marked_faces.count( face( opposite(candidate, tmesh), tmesh) ) ); - - uf.unify_sets( handles[hd], handles[opposite(candidate, tmesh)] ); - } - - std::size_t nb_cc = uf.number_of_sets(); - if ( nb_cc != 1 ) - { - // if more than one connected component is found then the patch - // made of marked faces contains "non-manifold" vertices. - // The smallest components need to be marked so that the patch - // made of marked faces is a topological disk - - // we will explore in parallel the connected components and will stop - // when all but one connected component have been entirely explored. - // We add one face at a time for each cc in order to not explore a - // potentially very large cc. - std::vector< std::vector > stacks_per_cc(nb_cc); - std::vector< std::set > faces_per_cc(nb_cc); - std::vector< bool > exploration_finished(nb_cc, false); - - // init the stacks of halfedges using the cc of the boundary - std::size_t index=0; - std::map< halfedge_descriptor, std::size_t > ccs; - typedef std::pair Pair_type; - for(Pair_type p : handles) - { - halfedge_descriptor opp_hedge = opposite(p.first, tmesh); - if (is_border(opp_hedge, tmesh)) continue; // nothing to do on the boundary - - typedef typename std::map< halfedge_descriptor, std::size_t >::iterator Map_it; - std::pair insert_res= - ccs.insert( std::make_pair(*uf.find( p.second ), index) ); - - if (insert_res.second) - ++index; - - stacks_per_cc[ insert_res.first->second ].push_back( prev(opp_hedge, tmesh) ); - stacks_per_cc[ insert_res.first->second ].push_back( next(opp_hedge, tmesh) ); - faces_per_cc[ insert_res.first->second ].insert( face(opp_hedge, tmesh) ); - } - - if (index != nb_cc) - { - // most probably, one cc is a cycle of border edges -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Trying to remove a component with a cycle of halfedges (nested hole or whole component), not handled yet.\n"; -#endif - all_removed=false; - continue; - } - - std::size_t nb_ccs_to_be_explored = nb_cc; - index=0; - //explore the cc's - do - { - // try to extract one more face for a given cc - do - { - CGAL_assertion( !exploration_finished[index] ); - CGAL_assertion( !stacks_per_cc.empty() ); - CGAL_assertion( !stacks_per_cc[index].empty() ); - halfedge_descriptor hd = stacks_per_cc[index].back(); - stacks_per_cc[index].pop_back(); - hd = opposite(hd, tmesh); - if ( !is_border(hd,tmesh) && !marked_faces.count(face(hd, tmesh) ) ) - { - if ( faces_per_cc[index].insert( face(hd, tmesh) ).second ) - { - stacks_per_cc[index].push_back( next(hd, tmesh) ); - stacks_per_cc[index].push_back( prev(hd, tmesh) ); - break; - } - } - if (stacks_per_cc[index].empty()) break; - } - while(true); - - // the exploration of a cc is finished when its stack is empty - exploration_finished[index]=stacks_per_cc[index].empty(); - if ( exploration_finished[index] ) - --nb_ccs_to_be_explored; - if ( nb_ccs_to_be_explored==1 ) - break; - while ( exploration_finished[(++index)%nb_cc] ); - index=index%nb_cc; - } - while(true); - - /// \todo use the area criteria? this means maybe continue exploration of larger cc - // mark faces of completetly explored cc - for (index=0; index< nb_cc; ++index) - { - if( exploration_finished[index] ) - { - for(face_descriptor fd : faces_per_cc[index]) - marked_faces.insert(fd); - } - } - } - - // make sure the selection is a topological disk (otherwise we need another treatment) - if (!is_selection_a_topological_disk(marked_faces, tmesh)) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Trying to handle a non-topological disk, do nothing\n"; -#endif - all_removed=false; - continue; - } - - some_removed = true; - // collect simplices to be removed - std::set vertices_to_keep; - std::set halfedges_to_keep; - for(halfedge_descriptor hd : border) - { - if ( !marked_faces.count(face(opposite(hd, tmesh), tmesh)) ) - { - halfedges_to_keep.insert( hd ); - vertices_to_keep.insert( target(hd, tmesh) ); - } - } - - // backup next,prev relationships to set after patch removal - std::vector< std::pair > next_prev_halfedge_pairs; - halfedge_descriptor first_border_hd=*( halfedges_to_keep.begin() ); - halfedge_descriptor current_border_hd=first_border_hd; - do - { - halfedge_descriptor prev_border_hd=current_border_hd; - current_border_hd=next(current_border_hd, tmesh); - while( marked_faces.count( face( opposite(current_border_hd, tmesh), tmesh) ) ) - current_border_hd=next(opposite(current_border_hd, tmesh), tmesh); - next_prev_halfedge_pairs.push_back( std::make_pair(prev_border_hd, current_border_hd) ); - } - while(current_border_hd!=first_border_hd); - - // collect vertices and edges to remove and do remove faces - std::set edges_to_remove; - std::set vertices_to_remove; - for(face_descriptor fd : marked_faces) - { - halfedge_descriptor hd=halfedge(fd, tmesh); - for(int i=0; i<3; ++i) - { - if ( !halfedges_to_keep.count(hd) ) - edges_to_remove.insert( edge(hd, tmesh) ); - - if ( !vertices_to_keep.count(target(hd,tmesh)) ) - vertices_to_remove.insert( target(hd,tmesh) ); - - hd=next(hd, tmesh); - } - - remove_face(fd, tmesh); - face_set.erase(fd); - } - - // remove vertices - for(vertex_descriptor vd : vertices_to_remove) - remove_vertex(vd, tmesh); - - // remove edges - for(edge_descriptor ed : edges_to_remove) - { - degenerate_edges_to_remove.erase(ed); - remove_edge(ed, tmesh); - } - - // add a new face, set all border edges pointing to it - // and update halfedge vertex of patch boundary vertices - face_descriptor new_face = add_face(tmesh); - typedef std::pair Pair_type; - for(const Pair_type& p : next_prev_halfedge_pairs) - { - set_face(p.first, new_face, tmesh); - set_next(p.first, p.second, tmesh); - set_halfedge(target(p.first, tmesh), p.first, tmesh); - } - - set_halfedge(new_face, first_border_hd, tmesh); - // triangulate the new face and update the coordinate of the central vertex - halfedge_descriptor new_hd=Euler::add_center_vertex(first_border_hd, tmesh); - put(vpmap, target(new_hd, tmesh), pt); - - for(halfedge_descriptor hd : halfedges_around_target(new_hd, tmesh)) - { - if(is_degenerate_edge(edge(hd, tmesh), tmesh, np)) - degenerate_edges_to_remove.insert(edge(hd, tmesh)); - - if(face(hd, tmesh) != GT::null_face() && is_degenerate_triangle_face(face(hd, tmesh), tmesh)) - face_set.insert(face(hd, tmesh)); - } - - CGAL_assertion( is_valid_polygon_mesh(tmesh) ); - } - } - } - - return all_removed; -} - -template -bool remove_degenerate_edges(const EdgeRange& edge_range, - TriangleMesh& tmesh, - const CGAL_PMP_NP_CLASS& np) -{ - std::set::face_descriptor> face_set; - return remove_degenerate_edges(edge_range, tmesh, face_set, np); -} - -template -bool remove_degenerate_edges(TriangleMesh& tmesh, - const CGAL_PMP_NP_CLASS& np) -{ - std::set::face_descriptor> face_set; - return remove_degenerate_edges(edges(tmesh), tmesh, face_set, np); -} - -template -bool remove_degenerate_edges(const EdgeRange& edge_range, - TriangleMesh& tmesh) -{ - std::set::face_descriptor> face_set; - return remove_degenerate_edges(edge_range, tmesh, face_set, parameters::all_default()); -} - -template -bool remove_degenerate_edges(TriangleMesh& tmesh) -{ - std::set::face_descriptor> face_set; - return remove_degenerate_edges(edges(tmesh), tmesh, face_set, parameters::all_default()); -} - -// \ingroup PMP_repairing_grp -// removes the degenerate faces from a triangulated surface mesh. -// A face is considered degenerate if two of its vertices share the same location, -// or more generally if all its vertices are collinear. -// -// @pre `CGAL::is_triangle_mesh(tmesh)` -// -// @tparam TriangleMesh a model of `FaceListGraph` and `MutableFaceGraph` -// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" -// -// @param tmesh the triangulated surface mesh to be repaired -// @param np optional \ref pmp_namedparameters "Named Parameters" described below -// -// \cgalNamedParamsBegin -// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. -// The type of this map is model of `ReadWritePropertyMap`. -// If this parameter is omitted, an internal property map for -// `CGAL::vertex_point_t` must be available in `TriangleMesh` -// \cgalParamEnd -// \cgalParamBegin{geom_traits} a geometric traits class instance. -// The traits class must provide the nested type `Point_3`, -// and the nested functors: -// - `Compare_distance_3` to compute the distance between 2 points -// - `Collinear_3` to check whether 3 points are collinear -// - `Less_xyz_3` to compare lexicographically two points -/// - `Equal_3` to check whether 2 points are identical. -/// For each functor Foo, a function `Foo foo_object()` must be provided. -// \cgalParamEnd -// \cgalNamedParamsEnd -// -// \todo the function might not be able to remove all degenerate faces. -// We should probably do something with the return type. -// -/// \return `true` if all degenerate faces were successfully removed, and `false` otherwise. -template -bool remove_degenerate_faces(const FaceRange& face_range, - TriangleMesh& tmesh, - const NamedParameters& np) -{ - CGAL_assertion(CGAL::is_triangle_mesh(tmesh)); - - using parameters::get_parameter; - using parameters::choose_parameter; - - typedef TriangleMesh TM; - typedef typename boost::graph_traits GT; - typedef typename GT::edge_descriptor edge_descriptor; - typedef typename GT::halfedge_descriptor halfedge_descriptor; - typedef typename GT::face_descriptor face_descriptor; - typedef typename GT::vertex_descriptor vertex_descriptor; - - typedef typename GetVertexPointMap::type VertexPointMap; - VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_property_map(vertex_point, tmesh)); - typedef typename GetGeomTraits::type Traits; - Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits()); - - typedef typename boost::property_traits::value_type Point_3; - typedef typename boost::property_traits::reference Point_ref; - - std::set degenerate_face_set; - degenerate_faces(face_range, tmesh, std::inserter(degenerate_face_set, degenerate_face_set.begin()), np); - - const std::size_t faces_size = faces(tmesh).size(); - - if(degenerate_face_set.empty()) - return true; - - if(degenerate_face_set.size() == faces_size) - { - clear(tmesh); - return true; - } - - // Sanitize the face range by adding adjacent degenerate faces - const std::size_t range_size = face_range.size(); - bool is_range_full_mesh = (range_size == faces_size); - if(!is_range_full_mesh) - { - std::list faces_to_visit(degenerate_face_set.begin(), degenerate_face_set.end()); - - while(!faces_to_visit.empty()) - { - face_descriptor fd = faces_to_visit.front(); - faces_to_visit.pop_front(); - - for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) - { - for(halfedge_descriptor inc_hd : halfedges_around_target(hd, tmesh)) - { - face_descriptor adj_fd = face(inc_hd, tmesh); - if(adj_fd == GT::null_face() || adj_fd == fd) - continue; - - if(is_degenerate_triangle_face(adj_fd, tmesh)) - { - if(degenerate_face_set.insert(adj_fd).second) - { - // successful insertion means we did not know about this face before - faces_to_visit.push_back(adj_fd); - } - } - } - } - } - } - - // Note that there can't be any null edge incident to the degenerate faces range, - // otherwise we would have a null face incident to the face range, and that is not possible - // after the sanitization above - std::set edge_range; - for(face_descriptor fd : degenerate_face_set) - for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) - edge_range.insert(edge(hd, tmesh)); - - // First remove edges of length 0 - bool all_removed = remove_degenerate_edges(edge_range, tmesh, degenerate_face_set, np); - - CGAL_assertion_code(for(face_descriptor fd : degenerate_face_set) {) - CGAL_assertion(is_degenerate_triangle_face(fd, tmesh)); - CGAL_assertion_code(}) - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - { - std::cout <<"Done with null edges.\n"; - std::ofstream output("/tmp/no_null_edges.off"); - output << std::setprecision(17) << tmesh << "\n"; - output.close(); - } -#endif - - // Then, remove triangles made of 3 collinear points - - // start by filtering out border faces - // TODO: shall we avoid doing that in case a non-manifold vertex on the boundary or if a whole component disappear? - std::set border_deg_faces; - for(face_descriptor f : degenerate_face_set) - { - halfedge_descriptor h = halfedge(f, tmesh); - for (int i=0; i<3; ++i) - { - if ( is_border( opposite(h, tmesh), tmesh) ) - { - border_deg_faces.insert(f); - break; - } - h = next(h, tmesh); - } - } - - while( !border_deg_faces.empty() ) - { - face_descriptor f_to_rm = *border_deg_faces.begin(); - border_deg_faces.erase(border_deg_faces.begin()); - - halfedge_descriptor h = halfedge(f_to_rm, tmesh); - for (int i=0; i<3; ++i) - { - face_descriptor f = face(opposite(h, tmesh), tmesh); - if ( f!=GT::null_face() ) - { - if (is_degenerate_triangle_face(f, tmesh, np) ) - border_deg_faces.insert(f); - } - h = next(h, tmesh); - } - - while( !is_border(opposite(h, tmesh), tmesh) ) - { - h = next(h, tmesh); - } - - degenerate_face_set.erase(f_to_rm); - Euler::remove_face(h, tmesh); - } - - // Ignore faces with null edges - if(!all_removed) - { - std::map are_degenerate_edges; - - for(face_descriptor fd : degenerate_face_set) - { - for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) - { - edge_descriptor ed = edge(hd, tmesh); - std::pair::iterator, bool> is_insert_successful = - are_degenerate_edges.insert(std::make_pair(ed, false)); - - bool is_degenerate = false; - if(is_insert_successful.second) - { - // did not previously exist in the map, so actually have to check if it is degenerate - if(traits.equal_3_object()(get(vpmap, target(ed, tmesh)), get(vpmap, source(ed, tmesh)))) - is_degenerate = true; - } - - is_insert_successful.first->second = is_degenerate; - - if(is_degenerate) - { - halfedge_descriptor h = halfedge(ed, tmesh); - if(!is_border(h, tmesh)) - degenerate_face_set.erase(face(h, tmesh)); - - h = opposite(h, tmesh); - if(!is_border(h, tmesh)) - degenerate_face_set.erase(face(h, tmesh)); - } - } - } - } - - // first remove degree 3 vertices that are part of a cap - // (only the vertex in the middle of the opposite edge) - // This removal does not change the shape of the mesh. - while (!degenerate_face_set.empty()) - { - std::set vertices_to_remove; - for(face_descriptor fd : degenerate_face_set) - { - for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh)) - { - vertex_descriptor vd = target(hd, tmesh); - if (degree(vd, tmesh) == 3 && is_border(vd, tmesh)==GT::null_halfedge()) - { - vertices_to_remove.insert(vd); - break; - } - } - } - - for(vertex_descriptor vd : vertices_to_remove) - { - halfedge_descriptor hd=halfedge(vd, tmesh); - for(halfedge_descriptor hd2 : halfedges_around_target(hd, tmesh)) - if (!is_border(hd2, tmesh)) - degenerate_face_set.erase( face(hd2, tmesh) ); - - // remove the central vertex and check if the new face is degenerated - hd=CGAL::Euler::remove_center_vertex(hd, tmesh); - if (is_degenerate_triangle_face(face(hd, tmesh), tmesh, np)) - { - degenerate_face_set.insert( face(hd, tmesh) ); - } - } - - if (vertices_to_remove.empty()) - break; - } - - while (!degenerate_face_set.empty()) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << "Loop on removing deg faces\n"; - // ensure the mesh is not broken - { - std::ofstream out("/tmp/out.off"); - out << tmesh; - out.close(); - - std::vector points; - std::vector > triangles; - std::ifstream in("/tmp/out.off"); - CGAL::read_OFF(in, points, triangles); - if (!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)) - { - std::cerr << "Warning: got a polygon soup (may simply be a non-manifold vertex)!\n"; - } - } -#endif - - face_descriptor fd = *degenerate_face_set.begin(); - - // look whether an incident triangle is also degenerate - bool detect_cc_of_degenerate_triangles = false; - for(halfedge_descriptor hd : halfedges_around_face(halfedge(fd, tmesh), tmesh) ) - { - face_descriptor adjacent_face = face( opposite(hd, tmesh), tmesh ); - if ( adjacent_face!=GT::null_face() && - degenerate_face_set.count(adjacent_face) ) - { - detect_cc_of_degenerate_triangles = true; - break; - } - } - - if (!detect_cc_of_degenerate_triangles) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " no degenerate neighbors, using a flip.\n"; -#endif - degenerate_face_set.erase(degenerate_face_set.begin()); - - // flip the longest edge of the triangle - Point_ref p1 = get(vpmap, target( halfedge(fd, tmesh), tmesh) ); - Point_ref p2 = get(vpmap, target(next(halfedge(fd, tmesh), tmesh), tmesh) ); - Point_ref p3 = get(vpmap, source( halfedge(fd, tmesh), tmesh) ); - - CGAL_assertion(p1!=p2 && p1!=p3 && p2!=p3); - - typename Traits::Compare_distance_3 compare_distance = traits.compare_distance_3_object(); - - halfedge_descriptor edge_to_flip; - if (compare_distance(p1,p2, p1,p3) != CGAL::SMALLER) // p1p2 > p1p3 - { - if (compare_distance(p1,p2, p2,p3) != CGAL::SMALLER) // p1p2 > p2p3 - // flip p1p2 - edge_to_flip = next( halfedge(fd, tmesh), tmesh ); - else - // flip p2p3 - edge_to_flip = prev( halfedge(fd, tmesh), tmesh ); - } - else - if (compare_distance(p1,p3, p2,p3) != CGAL::SMALLER) // p1p3>p2p3 - //flip p3p1 - edge_to_flip = halfedge(fd, tmesh); - else - //flip p2p3 - edge_to_flip = prev( halfedge(fd, tmesh), tmesh ); - - face_descriptor opposite_face=face( opposite(edge_to_flip, tmesh), tmesh); - if ( opposite_face == GT::null_face() ) - { - // simply remove the face - Euler::remove_face(edge_to_flip, tmesh); - } - else - { - // condition for the flip to be valid (the edge to be created do not already exists) - if ( !halfedge(target(next(edge_to_flip, tmesh), tmesh), - target(next(opposite(edge_to_flip, tmesh), tmesh), tmesh), - tmesh).second ) - { - Euler::flip_edge(edge_to_flip, tmesh); - } - else - { - all_removed=false; -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " WARNING: flip is not possible\n"; - // \todo Let p and q be the vertices opposite to `edge_to_flip`, and let - // r be the vertex of `edge_to_flip` that is the furthest away from - // the edge `pq`. In that case I think we should remove all the triangles - // so that the triangle pqr is in the mesh. -#endif - } - } - } - else - { - // Process a connected component of degenerate faces - // get all the faces from the connected component - // and the boundary edges - std::set cc_faces; - std::vector queue; - std::vector boundary_hedges; - std::vector inside_hedges; - queue.push_back(fd); - cc_faces.insert(fd); - - while(!queue.empty()) - { - face_descriptor top=queue.back(); - queue.pop_back(); - for(halfedge_descriptor hd : halfedges_around_face(halfedge(top, tmesh), tmesh) ) - { - face_descriptor adjacent_face = face( opposite(hd, tmesh), tmesh ); - if ( adjacent_face==GT::null_face() || - degenerate_face_set.count(adjacent_face)==0 ) - { - boundary_hedges.push_back(hd); - } - else - { - if (cc_faces.insert(adjacent_face).second) - queue.push_back(adjacent_face); - - if ( hd < opposite(hd, tmesh) ) - inside_hedges.push_back(hd); - } - } - } - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " Deal with a cc of " << cc_faces.size() << " degenerate faces.\n"; - /// dump cc_faces - { - int id=0; - std::map vids; - for(face_descriptor f : cc_faces) - { - if ( vids.insert( std::make_pair( target(halfedge(f, tmesh), tmesh), id) ).second ) ++id; - if ( vids.insert( std::make_pair( target(next(halfedge(f, tmesh), tmesh), tmesh), id) ).second ) ++id; - if ( vids.insert( std::make_pair( target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh), id) ).second ) ++id; - } - - std::ofstream output("/tmp/cc_faces.off"); - output << std::setprecision(17); - output << "OFF\n" << vids.size() << " " << cc_faces.size() << " 0\n"; - std::vector points(vids.size()); - typedef std::pair Pair_type; - for(Pair_type p : vids) - - points[p.second]=get(vpmap, p.first); - for(typename Traits::Point_3 p : points) - output << p << "\n"; - - for(face_descriptor f : cc_faces) - { - output << "3 " - << vids[ target(halfedge(f, tmesh), tmesh) ] << " " - << vids[ target(next(halfedge(f, tmesh), tmesh), tmesh) ] << " " - << vids[ target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh) ] << "\n"; - } - - for (std::size_t pid=2; pid!=points.size(); ++pid) - { - CGAL_assertion(collinear(points[0], points[1], points[pid])); - } - } -#endif - - // find vertices strictly inside the cc - std::set boundary_vertices; - for(halfedge_descriptor hd : boundary_hedges) - boundary_vertices.insert( target(hd, tmesh) ); - - std::set inside_vertices; - for(halfedge_descriptor hd : inside_hedges) - { - if (!boundary_vertices.count( target(hd, tmesh) )) - inside_vertices.insert( target(hd, tmesh) ); - if (!boundary_vertices.count( source(hd, tmesh) )) - inside_vertices.insert( source(hd, tmesh) ); - } - - // v-e+f = 1 for a topological disk and e = (3f+#boundary_edges)/2 - if (boundary_vertices.size()+inside_vertices.size() - - (cc_faces.size()+boundary_hedges.size())/2 != 1) - { - //cc_faces does not define a topological disk - /// \todo Find to way to handle that case -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " WARNING: Cannot remove the component of degenerate faces: not a topological disk.\n"; -#endif - - for(face_descriptor f : cc_faces) - degenerate_face_set.erase(f); - - continue; - } - - // preliminary step to check if the operation is possible - // sort the boundary points along the common supporting line - // we first need a reference point - typedef internal::Less_vertex_point Less_vertex; - std::pair< - typename std::set::iterator, - typename std::set::iterator > ref_vertices = - boost::minmax_element( boundary_vertices.begin(), - boundary_vertices.end(), - Less_vertex(traits, vpmap) ); - - // and then we sort the vertices using this reference point - typedef std::set Sorted_point_set; - Sorted_point_set sorted_points; - for(vertex_descriptor v : boundary_vertices) - sorted_points.insert( get(vpmap,v) ); - - CGAL_assertion(sorted_points.size()== - std::set(sorted_points.begin(), - sorted_points.end()).size()); - - CGAL_assertion( get( vpmap, *ref_vertices.first)==*sorted_points.begin() ); - CGAL_assertion( get( vpmap, *ref_vertices.second)==*std::prev(sorted_points.end()) ); - - const typename Traits::Point_3& xtrm1 = *sorted_points.begin(); - const typename Traits::Point_3& xtrm2 = *std::prev(sorted_points.end()); - - // recover halfedges on the hole, bounded by the reference vertices - std::vector side_one, side_two; - - // look for the outgoing border halfedge of the first extreme point - for(halfedge_descriptor hd : boundary_hedges) - { - if ( get(vpmap, source(hd, tmesh)) == xtrm1 ) - { - side_one.push_back(hd); - break; - } - } - CGAL_assertion(side_one.size()==1); - - bool non_monotone_border = false; - - while( get(vpmap, target(side_one.back(), tmesh)) != xtrm2 ) - { - vertex_descriptor prev_vertex = target(side_one.back(), tmesh); - for(halfedge_descriptor hd : boundary_hedges) - { - if ( source(hd, tmesh) == prev_vertex ) - { - if ( get(vpmap, target(hd, tmesh)) < get(vpmap, prev_vertex) ) - non_monotone_border = true; - - side_one.push_back(hd); - break; - } - } - - if (non_monotone_border) - break; - } - - if (non_monotone_border) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " WARNING: Cannot remove the component of degenerate faces: border not a monotonic cycle.\n"; -#endif - - for(face_descriptor f : cc_faces) - degenerate_face_set.erase(f); - - continue; - } - - // look for the outgoing border halfedge of second extreme vertex - for(halfedge_descriptor hd : boundary_hedges) - { - if ( source(hd, tmesh) == target(side_one.back(), tmesh) ) - { - side_two.push_back(hd); - break; - } - } - CGAL_assertion(side_two.size()==1); - - while( target(side_two.back(), tmesh) != source(side_one.front(), tmesh) ) - { - vertex_descriptor prev_vertex = target(side_two.back(), tmesh); - for(halfedge_descriptor hd : boundary_hedges) - { - if ( source(hd, tmesh) == prev_vertex ) - { - if ( get(vpmap, target(hd, tmesh)) > get(vpmap, prev_vertex) ) - non_monotone_border = true; - - side_two.push_back(hd); - break; - } - } - - if (non_monotone_border) - break; - } - - if (non_monotone_border) - { -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " WARNING: Cannot remove the component of degenerate faces: border not a monotonic cycle.\n"; -#endif - - for(face_descriptor f : cc_faces) - degenerate_face_set.erase(f); - - continue; - } - - CGAL_assertion( side_one.size()+side_two.size()==boundary_hedges.size() ); - - // reverse the order of the second side so as to follow - // the same order than side one - std::reverse(side_two.begin(), side_two.end()); - for(halfedge_descriptor& h : side_two) - h=opposite(h, tmesh); - - //make sure the points of the vertices along side_one are correctly sorted - std::vector side_points; - side_points.reserve(side_one.size()+1); - side_points.push_back(get(vpmap,source(side_one.front(), tmesh))); - - for(halfedge_descriptor h : side_one) - side_points.push_back(get(vpmap,target(h, tmesh))); - - CGAL_assertion(get(vpmap,source(side_one.front(), tmesh))==side_points.front()); - CGAL_assertion(get(vpmap,target(side_one.back(), tmesh))==side_points.back()); - - //\todo the reordering could lead to the apparition of null edges. - std::sort(side_points.begin(), side_points.end()); - - CGAL_assertion(std::unique(side_points.begin(), side_points.end())==side_points.end()); - for(std::size_t i=0;i::iterator side_one_it = side_one.begin(); - typename std::vector::iterator side_two_it = side_two.begin(); - for(;it_pt!=it_pt_end;++it_pt) - { - // check if it_pt is the point of the target of one or two halfedges - bool target_of_side_one = get(vpmap, target(*side_one_it, tmesh))==*it_pt; - bool target_of_side_two = get(vpmap, target(*side_two_it, tmesh))==*it_pt; - - if (target_of_side_one && target_of_side_two) - { - for(halfedge_descriptor h : halfedges_around_target(*side_one_it, tmesh)) - { - if (source(h, tmesh)==target(*side_two_it, tmesh)) - { - non_collapsable=true; - break; - } - } - } - else - { - CGAL_assertion(target_of_side_one || target_of_side_two); - vertex_descriptor v1 = target_of_side_one ? target(*side_one_it, tmesh) - : target(*side_two_it, tmesh); - vertex_descriptor v2 = target_of_side_two ? target(next(opposite(*side_one_it, tmesh), tmesh), tmesh) - : target(next(*side_two_it, tmesh), tmesh); - for(halfedge_descriptor h : halfedges_around_target(v1, tmesh)) - { - if (source(h, tmesh)==v2) - { - non_collapsable=true; - break; - } - } - } - - if(non_collapsable) break; - if (target_of_side_one) ++side_one_it; - if (target_of_side_two) ++side_two_it; - } - - if (non_collapsable) - { - for(face_descriptor f : cc_faces) - degenerate_face_set.erase(f); - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " WARNING: cannot remove a connected components of degenerate faces.\n"; -#endif - continue; - } - - // now proceed to the fix - // update the face and halfedge vertex pointers on the boundary - for(halfedge_descriptor h : boundary_hedges) - { - set_face(h, GT::null_face(), tmesh); - set_halfedge(target(h,tmesh), h, tmesh); - } - - // update next/prev pointers of boundary_hedges - for(halfedge_descriptor h : boundary_hedges) - { - halfedge_descriptor next_candidate = next( h, tmesh); - while (face(next_candidate, tmesh)!=GT::null_face()) - next_candidate = next( opposite( next_candidate, tmesh), tmesh); - - set_next(h, next_candidate, tmesh); - } - - // remove degenerate faces - for(face_descriptor f : cc_faces) - { - degenerate_face_set.erase(f); - remove_face(f, tmesh); - } - - // remove interior edges - for(halfedge_descriptor h : inside_hedges) - remove_edge(edge(h, tmesh), tmesh); - - // remove interior vertices - for(vertex_descriptor v : inside_vertices) - remove_vertex(v, tmesh); - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - std::cout << " side_one.size() " << side_one.size() << "\n"; - std::cout << " side_two.size() " << side_two.size() << "\n"; -#endif - - CGAL_assertion( source(side_one.front(), tmesh) == *ref_vertices.first ); - CGAL_assertion( source(side_two.front(), tmesh) == *ref_vertices.first ); - CGAL_assertion( target(side_one.back(), tmesh) == *ref_vertices.second ); - CGAL_assertion( target(side_two.back(), tmesh) == *ref_vertices.second ); - - // now split each side to contains the same sequence of points - // first side - int hi=0; - for (typename Sorted_point_set::iterator it=std::next(sorted_points.begin()), - it_end=sorted_points.end(); it!=it_end; ++it) - { - CGAL_assertion( *std::prev(it) == get(vpmap, source(side_one[hi], tmesh) ) ); - if( *it != get(vpmap, target(side_one[hi], tmesh) ) ) - { - // split the edge and update the point - halfedge_descriptor h1 = next(opposite(side_one[hi], tmesh), tmesh); - put(vpmap, - target(Euler::split_edge(side_one[hi], tmesh), tmesh), - *it); - - // split_edge updates the halfedge of the source vertex of h, - // since we reuse later the halfedge of the first refernce vertex - // we must set it as we need. - if ( source(h1,tmesh) == *ref_vertices.first) - set_halfedge(*ref_vertices.first, prev( prev(side_one[hi], tmesh), tmesh), tmesh ); - - // retriangulate the opposite face - if ( face(h1, tmesh) != GT::null_face()) - Euler::split_face(h1, opposite(side_one[hi], tmesh), tmesh); - } - else - { - ++hi; - } - } - - // second side - hi=0; - for (typename Sorted_point_set::iterator it=std::next(sorted_points.begin()), - it_end=sorted_points.end(); it!=it_end; ++it) - { - CGAL_assertion( *std::prev(it) == get(vpmap, source(side_two[hi], tmesh) ) ); - if( *it != get(vpmap, target(side_two[hi], tmesh) ) ) - { - // split the edge and update the point - halfedge_descriptor h2 = Euler::split_edge(side_two[hi], tmesh); - put(vpmap, target(h2, tmesh), *it); - - // split_edge updates the halfedge of the source vertex of h, - // since we reuse later the halfedge of the first refernce vertex - // we must set it as we need. - if ( source(h2,tmesh) == *ref_vertices.first) - set_halfedge(*ref_vertices.first, opposite( h2, tmesh), tmesh ); - - // retriangulate the face - if ( face(h2, tmesh) != GT::null_face()) - Euler::split_face(h2, next(side_two[hi], tmesh), tmesh); - } - else - { - ++hi; - } - } - - CGAL_assertion( target(halfedge(*ref_vertices.first, tmesh), tmesh) == *ref_vertices.first ); - CGAL_assertion( face(halfedge(*ref_vertices.first, tmesh), tmesh) == GT::null_face() ); - -#ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG - { - halfedge_descriptor h_side2 = halfedge(*ref_vertices.first, tmesh); - halfedge_descriptor h_side1 = next(h_side2, tmesh); - - do - { - CGAL_assertion( get(vpmap, source(h_side1, tmesh)) == get(vpmap, target(h_side2, tmesh)) ); - CGAL_assertion( get(vpmap, target(h_side1, tmesh)) == get(vpmap, source(h_side2, tmesh)) ); - - if ( target(next(opposite(h_side1, tmesh), tmesh), tmesh) == - target(next(opposite(h_side2, tmesh), tmesh), tmesh) ) - { - CGAL_error_msg("Forbidden simplification"); - } - - h_side2 = prev(h_side2, tmesh); - h_side1 = next(h_side1, tmesh); - } - while( target(h_side1, tmesh) != *ref_vertices.second ); - } -#endif - - // remove side1 and replace its opposite hedges by those of side2 - halfedge_descriptor h_side2 = halfedge(*ref_vertices.first, tmesh); - halfedge_descriptor h_side1 = next(h_side2, tmesh); - while(true) - { - CGAL_assertion( get(vpmap, source(h_side1, tmesh)) == get(vpmap, target(h_side2, tmesh)) ); - CGAL_assertion( get(vpmap, target(h_side1, tmesh)) == get(vpmap, source(h_side2, tmesh)) ); - - // backup target vertex - vertex_descriptor vertex_to_remove = target(h_side1, tmesh); - if (vertex_to_remove!=*ref_vertices.second) - { - vertex_descriptor replacement_vertex = source(h_side2, tmesh); - // replace the incident vertex - for(halfedge_descriptor hd : halfedges_around_target(h_side1, tmesh)) - set_target(hd, replacement_vertex, tmesh); - } - - // prev side2 hedge for next loop - halfedge_descriptor h_side2_for_next_turn = prev(h_side2, tmesh); - - // replace the opposite of h_side1 by h_side2 - halfedge_descriptor opposite_h_side1 = opposite( h_side1, tmesh); - face_descriptor the_face = face(opposite_h_side1, tmesh); - set_face(h_side2, the_face, tmesh); - - if (the_face!=GT::null_face()) - set_halfedge(the_face, h_side2, tmesh); - - set_next(h_side2, next(opposite_h_side1, tmesh), tmesh); - set_next(prev(opposite_h_side1, tmesh), h_side2, tmesh); - - // take the next hedges - edge_descriptor edge_to_remove = edge(h_side1, tmesh); - h_side1 = next(h_side1, tmesh); - - // now remove the extra edge - remove_edge(edge_to_remove, tmesh); - - // ... and the extra vertex if it's not the second reference - if (vertex_to_remove==*ref_vertices.second) - { - // update the halfedge pointer of the last vertex (others were already from side 2) - CGAL_assertion( target(opposite(h_side2, tmesh), tmesh) == vertex_to_remove ); - set_halfedge(vertex_to_remove, opposite(h_side2, tmesh), tmesh); - break; - } - else - { - remove_vertex(vertex_to_remove , tmesh); - } - - h_side2 = h_side2_for_next_turn; - } - } - } - - return all_removed; -} - -template -bool remove_degenerate_faces(const FaceRange& face_range, - TriangleMesh& tmesh) -{ - return remove_degenerate_faces(face_range, tmesh, parameters::all_default()); -} - -template -bool remove_degenerate_faces(TriangleMesh& tmesh, - const CGAL_PMP_NP_CLASS& np) -{ - return remove_degenerate_faces(faces(tmesh), tmesh, np); -} - -template -bool remove_degenerate_faces(TriangleMesh& tmesh) -{ - return remove_degenerate_faces(tmesh, - CGAL::Polygon_mesh_processing::parameters::all_default()); -} - -/// \ingroup PMP_repairing_grp -/// checks whether a vertex of a polygon mesh is non-manifold. -/// -/// @tparam PolygonMesh a model of `HalfedgeListGraph` -/// -/// @param v a vertex of `pm` -/// @param pm a triangle mesh containing `v` -/// -/// \warning This function has linear runtime with respect to the size of the mesh. -/// -/// \sa `duplicate_non_manifold_vertices()` -/// -/// \return `true` if the vertex is non-manifold, `false` otherwise. -template -bool is_non_manifold_vertex(typename boost::graph_traits::vertex_descriptor v, - const PolygonMesh& pm) -{ - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - typedef CGAL::dynamic_halfedge_property_t Halfedge_property_tag; - typedef typename boost::property_map::const_type Visited_halfedge_map; - - // Dynamic pmaps do not have default initialization values (yet) - Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm); - for(halfedge_descriptor h : halfedges(pm)) - put(visited_halfedges, h, false); - - std::size_t incident_null_faces_counter = 0; - for(halfedge_descriptor h : halfedges_around_target(v, pm)) - { - put(visited_halfedges, h, true); - if(CGAL::is_border(h, pm)) - ++incident_null_faces_counter; - } - - if(incident_null_faces_counter > 1) - { - // The vertex is the sole connection between two connected components --> non-manifold - return true; - } - - for(halfedge_descriptor h : halfedges(pm)) - { - if(v == target(h, pm)) - { - // Haven't seen that halfedge yet ==> more than one umbrella incident to 'v' ==> non-manifold - if(!get(visited_halfedges, h)) - return true; - } - } - - return false; -} - -namespace internal { - -template -struct Vertex_collector -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - bool has_old_vertex(const vertex_descriptor v) const { return collections.count(v) != 0; } - void tag_old_vertex(const vertex_descriptor v) - { - CGAL_precondition(!has_old_vertex(v)); - collections[v]; - } - - void collect_vertices(vertex_descriptor v1, vertex_descriptor v2) - { - std::vector& verts = collections[v1]; - if(verts.empty()) - verts.push_back(v1); - verts.push_back(v2); - } - - template - void dump(OutputIterator out) - { - typedef std::pair > Pair_type; - for(const Pair_type& p : collections) - *out++ = p.second; - } - - void dump(Emptyset_iterator) { } - - std::map > collections; -}; - -template -typename boost::graph_traits::vertex_descriptor -create_new_vertex_for_sector(typename boost::graph_traits::halfedge_descriptor sector_begin_h, - typename boost::graph_traits::halfedge_descriptor sector_last_h, - PolygonMesh& pm, - const VPM& vpm, - const ConstraintMap& cmap) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - vertex_descriptor old_vd = target(sector_begin_h, pm); - vertex_descriptor new_vd = add_vertex(pm); - put(vpm, new_vd, get(vpm, old_vd)); - - put(cmap, new_vd, true); - - set_halfedge(new_vd, sector_begin_h, pm); - halfedge_descriptor h = sector_begin_h; - do - { - set_target(h, new_vd, pm); - - if(h == sector_last_h) - break; - else - h = prev(opposite(h, pm), pm); - } - while(h != sector_begin_h); // for safety - CGAL_assertion(h != sector_begin_h); - - return new_vd; -} - -template -std::size_t make_umbrella_manifold(typename boost::graph_traits::halfedge_descriptor h, - PolygonMesh& pm, - internal::Vertex_collector& dmap, - const NamedParameters& np) -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - using parameters::get_parameter; - using parameters::choose_parameter; - - typedef typename GetVertexPointMap::type VertexPointMap; - VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_property_map(vertex_point, pm)); - - typedef typename internal_np::Lookup_named_param_def // default (no constraint pmap) - >::type VerticesMap; - VerticesMap cmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), - Constant_property_map(false)); - - std::size_t nb_new_vertices = 0; - - vertex_descriptor old_v = target(h, pm); - put(cmap, old_v, true); // store the duplicates - - // count the number of borders - int border_counter = 0; - halfedge_descriptor ih = h, done = ih, border_h = h; - do - { - if(is_border(ih, pm)) - { - border_h = ih; - ++border_counter; - } - - ih = prev(opposite(ih, pm), pm); - } - while(ih != done); - - const bool is_non_manifold_within_umbrella = (border_counter > 1); - if(!is_non_manifold_within_umbrella) - { - const bool first_time_meeting_v = !dmap.has_old_vertex(old_v); - if(first_time_meeting_v) - { - // The star is manifold, so if it is the first time we have met that vertex, - // there is nothing to do, we just keep the same vertex. - set_halfedge(old_v, h, pm); // to ensure halfedge(old_v, pm) stays valid - dmap.tag_old_vertex(old_v); // so that we know we have met old_v already, next time, we'll have to duplicate - } - else - { - // This is not the canonical star associated to 'v'. - // Create a new vertex, and move the whole star to that new vertex - halfedge_descriptor last_h = opposite(next(h, pm), pm); - vertex_descriptor new_v = create_new_vertex_for_sector(h, last_h, pm, vpm, cmap); - dmap.collect_vertices(old_v, new_v); - nb_new_vertices = 1; - } - } - // if there is more than one sector, look at each sector and split them away from the main one - else - { - // the first manifold sector, described by two halfedges - halfedge_descriptor sector_start_h = border_h; - CGAL_assertion(is_border(border_h, pm)); - - bool should_stop = false; - bool is_main_sector = true; - do - { - CGAL_assertion(is_border(sector_start_h, pm)); - - // collect the sector and split it away if it must be - halfedge_descriptor sector_last_h = sector_start_h; - do - { - halfedge_descriptor next_h = prev(opposite(sector_last_h, pm), pm); - - if(is_border(next_h, pm)) - break; - - sector_last_h = next_h; - } - while(sector_last_h != sector_start_h); - CGAL_assertion(!is_border(sector_last_h, pm)); - CGAL_assertion(sector_last_h != sector_start_h); - - halfedge_descriptor next_start_h = prev(opposite(sector_last_h, pm), pm); - - // there are multiple CCs incident to this particular vertex, and we should create a new vertex - // if it's not the first umbrella around 'old_v' or not the first sector, but only not if it's - // both the first umbrella and first sector. - const bool must_create_new_vertex = (!is_main_sector || dmap.has_old_vertex(old_v)); - - // In any case, we must set up the next pointer correctly - set_next(sector_start_h, opposite(sector_last_h, pm), pm); - - if(must_create_new_vertex) - { - vertex_descriptor new_v = create_new_vertex_for_sector(sector_start_h, sector_last_h, pm, vpm, cmap); - dmap.collect_vertices(old_v, new_v); - ++nb_new_vertices; - } - else - { - // Ensure that halfedge(old_v, pm) stays valid - set_halfedge(old_v, sector_start_h, pm); - } - - is_main_sector = false; - sector_start_h = next_start_h; - should_stop = (sector_start_h == border_h); - } - while(!should_stop); - } - - return nb_new_vertices; -} - -} // end namespace internal - -/// \ingroup PMP_repairing_grp -/// collects the non-manifold vertices (if any) present in the mesh. A non-manifold vertex `v` is returned -/// via one incident halfedge `h` such that `target(h, pm) = v` for all the umbrellas that `v` apppears in -/// (an umbrella being the set of faces incident to all the halfedges reachable by walking around `v` -/// using `hnext = prev(opposite(h, pm), pm)`, starting from `h`). -/// -/// @tparam PolygonMesh a model of `HalfedgeListGraph` -/// @tparam OutputIterator a model of `OutputIterator` holding objects of type -/// `boost::graph_traits::%halfedge_descriptor` -/// -/// @param pm a triangle mesh -/// @param out the output iterator that collects halfedges incident to `v` -/// -/// \sa `is_non_manifold_vertex()` -/// \sa `duplicate_non_manifold_vertices()` -/// -/// \return the output iterator. -template -OutputIterator non_manifold_vertices(const PolygonMesh& pm, - OutputIterator out) -{ - // 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. - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - typedef CGAL::dynamic_vertex_property_t Vertex_bool_tag; - typedef typename boost::property_map::const_type Known_manifold_vertex_map; - typedef CGAL::dynamic_vertex_property_t Vertex_halfedge_tag; - typedef typename boost::property_map::const_type Visited_vertex_map; - typedef CGAL::dynamic_halfedge_property_t Halfedge_property_tag; - typedef typename boost::property_map::const_type Visited_halfedge_map; - - Known_manifold_vertex_map known_nm_vertices = get(Vertex_bool_tag(), pm); - Visited_vertex_map visited_vertices = get(Vertex_halfedge_tag(), pm); - Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm); - - const halfedge_descriptor null_h = boost::graph_traits::null_halfedge(); - - // Dynamic pmaps do not have default initialization values (yet) - for(vertex_descriptor v : vertices(pm)) - { - put(known_nm_vertices, v, false); - put(visited_vertices, v, null_h); - } - for(halfedge_descriptor h : halfedges(pm)) - put(visited_halfedges, h, false); - - for(halfedge_descriptor h : halfedges(pm)) - { - // If 'h' is not visited yet, we walk around the target of 'h' and mark these - // halfedges as visited. Thus, if we are here and the target is already marked as visited, - // it means that the vertex is non manifold. - if(!get(visited_halfedges, h)) - { - put(visited_halfedges, h, true); - bool is_non_manifold = false; - - vertex_descriptor v = target(h, pm); - if(get(visited_vertices, v) != null_h) // already seen this vertex, but not from this star - { - is_non_manifold = true; - // if this is the second time we visit that vertex and the first star was manifold, we have - // never reported the first star, but we must now - if(!get(known_nm_vertices, v)) - *out++ = get(visited_vertices, v); // that's a halfedge of the first star we've seen 'v' in - } - else - { - // first time we meet this vertex, just mark it so, and keep the halfedge we found the vertex with in memory - put(visited_vertices, v, h); - } - - // While walking the star of this halfedge, if we meet a border halfedge more than once, - // it means the mesh is pinched and we are also in the case of a non-manifold situation - halfedge_descriptor ih = h, done = ih; - int border_counter = 0; - do - { - put(visited_halfedges, ih, true); - if(is_border(ih, pm)) - ++border_counter; - - ih = prev(opposite(ih, pm), pm); - } - while(ih != done); - - if(border_counter > 1) - is_non_manifold = true; - - if(is_non_manifold) - { - *out++ = h; - put(known_nm_vertices, v, true); - } - } - } - - return out; -} - -/// \ingroup PMP_repairing_grp -/// duplicates all the non-manifold vertices of the input mesh. -/// -/// @tparam PolygonMesh a model of `HalfedgeListGraph` and `MutableHalfedgeGraph` -/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" -/// -/// @param pm the surface mesh to be repaired -/// @param np optional \ref pmp_namedparameters "Named Parameters" described below -/// -/// \cgalNamedParamsBegin -/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. -/// The type of this map is model of `ReadWritePropertyMap`. -/// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `PolygonMesh` -/// \cgalParamEnd -/// \cgalParamBegin{vertex_is_constrained_map} a writable property map with `vertex_descriptor` -/// as key and `bool` as `value_type`. `put(pmap, v, true)` will be called for each duplicated -/// vertices, as well as the original non-manifold vertex in the input mesh. -/// \cgalParamEnd -/// \cgalParamBegin{output_iterator} a model of `OutputIterator` with value type -/// `std::vector`. The first vertex of each vector is a non-manifold vertex -/// of the input mesh, followed by the new vertices that were created to fix this precise -/// non-manifold configuration. -/// \cgalParamEnd -/// \cgalNamedParamsEnd -/// -/// \return the number of vertices created. -template -std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm, - const NamedParameters& np) -{ - using parameters::get_parameter; - using parameters::choose_parameter; - - typedef boost::graph_traits GT; - typedef typename GT::halfedge_descriptor halfedge_descriptor; - - typedef typename internal_np::Lookup_named_param_def < - internal_np::output_iterator_t, - NamedParameters, - Emptyset_iterator - > ::type Output_iterator; - - Output_iterator out = choose_parameter(get_parameter(np, internal_np::output_iterator), Emptyset_iterator()); - - std::vector non_manifold_cones; - non_manifold_vertices(pm, std::back_inserter(non_manifold_cones)); - - internal::Vertex_collector dmap; - std::size_t nb_new_vertices = 0; - if(!non_manifold_cones.empty()) - { - for(halfedge_descriptor h : non_manifold_cones) - nb_new_vertices += internal::make_umbrella_manifold(h, pm, dmap, np); - - dmap.dump(out); - } - - return nb_new_vertices; -} - -template -std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm) -{ - return duplicate_non_manifold_vertices(pm, parameters::all_default()); -} - -/// \cond SKIP_IN_MANUAL -template -std::pair< bool, bool > -remove_self_intersections_one_step(TriangleMesh& tm, - std::set& faces_to_remove, - VertexPointMap& vpmap, - int step, - bool preserve_genus, - bool verbose) -{ - std::set faces_to_remove_copy = faces_to_remove; - - if (verbose) - std::cout << "DEBUG: running remove_self_intersections_one_step, step " << step - << " with " << faces_to_remove.size() << " intersecting faces\n"; - - CGAL_assertion(tm.is_valid()); - - typedef boost::graph_traits graph_traits; - typedef typename graph_traits::vertex_descriptor vertex_descriptor; - typedef typename graph_traits::edge_descriptor edge_descriptor; - typedef typename graph_traits::halfedge_descriptor halfedge_descriptor; - - bool something_was_done = false; // indicates if a region was successfully remeshed - bool all_fixed = true; // indicates if all removal went well - // indicates if a removal was not possible because the region handle has - // some boundary cycle of halfedges - bool topology_issue = false; - - if (verbose) - { - std::cout << " DEBUG: is_valid in one_step(tm)? "; - std::cout.flush(); - std::cout << is_valid_polygon_mesh(tm) << "\n"; - } - - if(!faces_to_remove.empty()){ - - while(!faces_to_remove.empty()) - { - // Process a connected component of faces to remove. - // collect all the faces from the connected component - std::set cc_faces; - std::vector queue(1, *faces_to_remove.begin()); // temporary queue - cc_faces.insert(queue.back()); - while(!queue.empty()) - { - face_descriptor top=queue.back(); - queue.pop_back(); - halfedge_descriptor h = halfedge(top,tm); - for (int i=0;i<3; ++i) - { - face_descriptor adjacent_face = face( opposite(h, tm), tm ); - if ( adjacent_face!=boost::graph_traits::null_face()) - { - if (faces_to_remove.count(adjacent_face) != 0 && - cc_faces.insert(adjacent_face).second) - queue.push_back(adjacent_face); - } - h = next(h, tm); - } - } - - // expand the region to be filled - if (step > 0) - expand_face_selection(cc_faces, tm, step, - make_boolean_property_map(cc_faces), - Emptyset_iterator()); - - // try to compactify the selection region by also selecting all the faces included - // in the bounding box of the initial selection - std::vector stack_for_expension; - Bbox_3 bb; - for(face_descriptor fd : cc_faces) - { - for(halfedge_descriptor h : halfedges_around_face(halfedge(fd, tm), tm)) - { - bb += get(vpmap, target(h, tm)).bbox(); - face_descriptor nf = face(opposite(h, tm), tm); - if (nf != boost::graph_traits::null_face() && - cc_faces.count(nf)==0) - { - stack_for_expension.push_back(opposite(h, tm)); - } - } - } - - while(!stack_for_expension.empty()) - { - halfedge_descriptor h=stack_for_expension.back(); - stack_for_expension.pop_back(); - if ( cc_faces.count(face(h,tm))==1) continue; - if ( do_overlap(bb, get(vpmap, target(next(h, tm), tm)).bbox()) ) - { - cc_faces.insert(face(h,tm)); - halfedge_descriptor candidate = opposite(next(h, tm), tm); - if ( face(candidate, tm) != boost::graph_traits::null_face() ) - stack_for_expension.push_back( candidate ); - candidate = opposite(prev(h, tm), tm); - if ( face(candidate, tm) != boost::graph_traits::null_face() ) - stack_for_expension.push_back( candidate ); - } - } - - // remove faces from the set to process - for(face_descriptor f : cc_faces) - faces_to_remove.erase(f); - - if (cc_faces.size()==1) continue; // it is a triangle nothing better can be done - - //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,tm); - for (int i=0;i<3; ++i) - { - if ( is_border( opposite(h, tm), tm) || - cc_faces.count( face( opposite(h, tm), tm) ) == 0) - { - boundary_hedges.push_back(h); - } - h=next(h, tm); - } - } - - // 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,tm)).second){ - bool any_face_added = false; - for(halfedge_descriptor hh : halfedges_around_target(h,tm)){ - if (!is_border(hh, tm)) - { - any_face_added |= cc_faces.insert(face(hh, tm)).second; // add the face to the current selection - faces_to_remove.erase(face(hh, tm)); - } - } - 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; - if(verbose) - std::cout << " DEBUG: CC not handled due to the presence at least one non-manifold vertex\n"; - 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, tm); - for (int i=0; i<3;++i) - { - if ( is_border(opposite(h, tm), tm) || - cc_faces.count( face(opposite(h, tm), tm) )== 0) - { - cc_border_hedges.push_back(h); - } - h=next(h, tm); - } - } - - if(!is_selection_a_topological_disk(cc_faces, tm)) - { - // check if the selection contains cycles of border halfedges - bool only_border_edges = true; - std::set mesh_border_hedge; - - for(halfedge_descriptor h : cc_border_hedges) - { - if ( !is_border(opposite(h, tm), tm) ) - only_border_edges = false; - else - mesh_border_hedge.insert( opposite(h, tm) ); - } - int nb_cycles=0; - 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, tm); - if (h==h_b) - { - // found a cycle - ++nb_cycles; - break; - } - else - { - typename std::set::iterator it = - mesh_border_hedge.find(h); - if ( it == mesh_border_hedge.end() ) - break; // not a cycle - mesh_border_hedge.erase(it); - } - }while(true); - } - - if(nb_cycles > (only_border_edges ? 1 : 0) ) - { - if(verbose) - std::cout << " DEBUG: CC not handled due to the presence of " - << nb_cycles << " of boundary edges\n"; - topology_issue = true; - continue; - } - else - { - if (preserve_genus) - { - if(verbose) - std::cout << " DEBUG: CC not handled because it is not a topological disk (preserve_genus=true)\n"; - all_fixed = false; - continue; - } - // count the number of cycles of halfedges of the boundary - std::map bhs; - for(halfedge_descriptor h : cc_border_hedges) - { - bhs[source(h, tm)]=target(h, tm); - } - 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){ - if(verbose) - std::cout << " DEBUG: CC not handled because it is not a topological disk(" - << nbc << " boundary cycles)\n"; - all_fixed = false; - continue; - } - else - { - if(verbose) - std::cout << " DEBUG: CC that is not a topological disk but has only one boundary cycle(preserve_genus=false)\n"; - } - } - } - - // sort halfedges so that they describe the sequence - // of halfedges of the hole to be made - CGAL_assertion( cc_border_hedges.size() > 2 ); - for(std::size_t i=0; i < cc_border_hedges.size()-2; ++i) - { - vertex_descriptor tgt = target(cc_border_hedges[i], tm); - for(std::size_t j=i+1; j cc_interior_vertices; - std::set cc_interior_edges; - - // first collect all vertices and edges incident to the faces to remove - for(face_descriptor fh : cc_faces) - { - for(halfedge_descriptor h : halfedges_around_face(halfedge(fh,tm),tm)) - { - if (halfedge(target(h, tm), tm)==h) // limit the number of insertions - cc_interior_vertices.insert(target(h, tm)); - cc_interior_edges.insert(edge(h,tm)); - } - } - // and then remove those on the boundary - for(halfedge_descriptor h : cc_border_hedges) - { - cc_interior_vertices.erase(target(h, tm)); - cc_interior_edges.erase(edge(h,tm)); - } - - if (verbose) - { - std::cout << " DEBUG: is_valid(tm) in one_step, before mesh changes? "; - std::cout << is_valid_polygon_mesh(tm) << std::endl; - } - - //try hole_filling. - typedef CGAL::Triple Face_indices; - typedef typename boost::property_traits::value_type Point; - std::vector hole_points, third_points; - hole_points.reserve(cc_border_hedges.size()); - third_points.reserve(cc_border_hedges.size()); - std::vector border_vertices; - for(halfedge_descriptor h : cc_border_hedges) - { - vertex_descriptor v = source(h, tm); - hole_points.push_back( get(vpmap, v) ); - border_vertices.push_back(v); - third_points.push_back(get(vpmap, target(next(opposite(h, tm), tm), tm))); // TODO fix me for mesh border edges - } - CGAL_assertion(hole_points.size() >= 3); - - // 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 patch; - if (hole_points.size()>3) - triangulate_hole_polyline(hole_points, - third_points, - std::back_inserter(patch)); - else - patch.push_back(Face_indices(0,1,2)); // trivial hole filling - - if(patch.empty()) - { -#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3 - if (verbose) - std::cout << " DEBUG: Failed to fill a hole using Delaunay search space.\n"; - triangulate_hole_polyline(hole_points, - third_points, - std::back_inserter(patch), - parameters::use_delaunay_triangulation(false)); -#endif - if (patch.empty()) - { - if (verbose) - std::cout << " DEBUG: Failed to fill a hole using the whole search space.\n"; - all_fixed = false; - continue; - } - } - - // make sure that the hole filling is valid, we check that no - // edge already in the mesh is present in patch. - bool non_manifold_edge_found = false; - for(const Face_indices& triangle : 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) - { - int vi=edges[2*k], vj=edges[2*k+1]; - // ignore boundary edges - if (vi+1==vj || (vj==0 && static_cast(vi)==border_vertices.size()-1) ) - continue; - halfedge_descriptor h = halfedge(border_vertices[vi], border_vertices[vj], tm).first; - if (h!=boost::graph_traits::null_halfedge() && - cc_interior_edges.count(edge(h, tm))==0) - { - non_manifold_edge_found=true; - break; - } - } - if (non_manifold_edge_found) break; - } - if (non_manifold_edge_found) - { - if (verbose) - std::cout << " DEBUG: Triangulation produced is non-manifold when plugged into the mesh.\n"; - all_fixed = false; - continue; - } - - // plug the new triangles in the mesh, reusing previous edges and faces - std::vector edge_stack(cc_interior_edges.begin(), cc_interior_edges.end()); - std::vector face_stack(cc_faces.begin(), cc_faces.end()); - - std::map< std::pair, halfedge_descriptor > halfedge_map; - int i=0; - // register border halfedges - for(halfedge_descriptor h : cc_border_hedges) - { - int j = static_cast( std::size_t(i+1)%cc_border_hedges.size() ); - halfedge_map.insert(std::make_pair( std::make_pair(i, j), h) ); - set_halfedge(target(h, tm), h, tm); // update vertex halfedge pointer - CGAL_assertion( border_vertices[i] == source(h, tm) && - border_vertices[j] == target(h, tm) ); - ++i; - } - - std::vector hedges; - hedges.reserve(4); - face_descriptor f = boost::graph_traits::null_face(); - for(const Face_indices& triangle : patch) - { - // get the new face - if (face_stack.empty()) - f=add_face(tm); - else - { - f=face_stack.back(); - face_stack.pop_back(); - } - - std::array indices = - make_array( triangle.first, - triangle.second, - triangle.third, - triangle.first ); - for (int i=0; i<3; ++i) - { - // get the corresponding halfedge (either a new one or an already created) - typename std::map< std::pair , halfedge_descriptor >::iterator insert_res = - halfedge_map.insert( - std::make_pair( std::make_pair(indices[i], indices[i+1]), - boost::graph_traits::null_halfedge() ) ).first; - if (insert_res->second == boost::graph_traits::null_halfedge()) - { - if (edge_stack.empty()) - insert_res->second = halfedge(add_edge(tm), tm); - else - { - insert_res->second = halfedge(edge_stack.back(), tm); - edge_stack.pop_back(); - } - - halfedge_map[std::make_pair(indices[i+1], indices[i])] = - opposite(insert_res->second, tm); - } - hedges.push_back(insert_res->second); - } - hedges.push_back(hedges.front()); - // update halfedge connections + face pointers - for(int i=0; i<3;++i) - { - set_next(hedges[i], hedges[i+1], tm); - set_face(hedges[i], f, tm); - set_target(hedges[i], border_vertices[indices[i+1]], tm); - } - set_halfedge(f, hedges[0], tm); - hedges.clear(); - } - - // now remove remaining edges, - for(edge_descriptor e : edge_stack) - remove_edge(e, tm); - // vertices, - for(vertex_descriptor vh : cc_interior_vertices) - remove_vertex(vh, tm); - // and remaning faces - for(face_descriptor f : face_stack) - remove_face(f, tm); - - if (verbose) - std::cout << " DEBUG: " << cc_faces.size() << " triangles removed, " - << patch.size() << " created\n"; - - CGAL_assertion(is_valid_polygon_mesh(tm)); - - something_was_done = true; - } - } - if (!something_was_done) - { - faces_to_remove.swap(faces_to_remove_copy); - if (verbose) - std::cout<<" DEBUG: Nothing was changed during this step, self-intersections won't be recomputed."< -bool remove_self_intersections(TriangleMesh& tm, const NamedParameters& np) -{ - typedef boost::graph_traits graph_traits; - typedef typename graph_traits::face_descriptor face_descriptor; - - // named parameter extraction - typedef typename GetVertexPointMap::type VertexPointMap; - VertexPointMap vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point), - get_property_map(vertex_point, tm)); - - const int max_steps = parameters::choose_parameter(parameters::get_parameter(np, internal_np::number_of_iterations), 7); - bool verbose = parameters::choose_parameter(parameters::get_parameter(np, internal_np::verbosity_level), 0) > 0; - bool preserve_genus = parameters::choose_parameter(parameters::get_parameter(np, internal_np::preserve_genus), true); - - if (verbose) - std::cout << "DEBUG: Starting remove_self_intersections, is_valid(tm)? " << is_valid_polygon_mesh(tm) << "\n"; - - // first handle the removal of degenerate faces - remove_degenerate_faces(tm, np); - - if (!preserve_genus) - duplicate_non_manifold_vertices(tm, np); - - if (verbose) - std::cout << "DEBUG: After degenerate faces removal, is_valid(tm)? " << is_valid_polygon_mesh(tm) << "\n"; - - // Look for self-intersections in the polyhedron and remove them - 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; - while( ++step Face_pair; - std::vector self_inter; - // TODO : possible optimization to reduce the range to check with the bbox - // of the previous patches or something. - self_intersections(tm, std::back_inserter(self_inter)); - - for(Face_pair fp : self_inter) - { - faces_to_remove.insert(fp.first); - faces_to_remove.insert(fp.second); - } - } - - if ( faces_to_remove.empty() && all_fixed){ - if (verbose) - std::cout<<"DEBUG: There is no more face to remove."< -bool remove_self_intersections(TriangleMesh& tm) -{ - return remove_self_intersections(tm, parameters::all_default()); -} -/// \endcond - -} } // end of CGAL::Polygon_mesh_processing +} // namespace Polygon_mesh_processing +} // namespace CGAL #endif // CGAL_POLYGON_MESH_PROCESSING_REPAIR_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h index 409bdc4900d..3272b946ff6 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h @@ -69,6 +69,7 @@ struct Intersect_facets OutputIterator* m_iterator; bool* m_intersected; }; + // typedefs typedef typename Kernel::Segment_3 Segment; typedef typename Kernel::Triangle_3 Triangle; @@ -76,7 +77,7 @@ struct Intersect_facets typedef typename boost::property_map::const_type Ppmap; // members - const TM& m_tmesh; + const TM& g; const VertexPointMap m_vpmap; mutable OutputIterator m_iterator; mutable bool m_intersected; @@ -86,10 +87,9 @@ struct Intersect_facets typename Kernel::Construct_triangle_3 triangle_functor; typename Kernel::Do_intersect_3 do_intersect_3_functor; - Intersect_facets(const TM& tmesh, OutputIterator it, VertexPointMap vpmap, const Kernel& kernel) : - m_tmesh(tmesh), + g(tmesh), m_vpmap(vpmap), m_iterator(it), m_intersected(false), @@ -99,100 +99,115 @@ struct Intersect_facets do_intersect_3_functor(kernel.do_intersect_3_object()) { } + halfedge_descriptor common_vertex(halfedge_descriptor& bh, halfedge_descriptor ch) const + { + halfedge_descriptor v = boost::graph_traits::null_halfedge(); + + if(target(bh, g) == target(ch, g)) + v = ch; + else if(target(bh, g) == target(next(ch, g), g)) + v = next(ch, g); + else if(target(bh, g) == target(next(next(ch, g), g), g)) + v = next(next(ch, g), g); + + if(v == boost::graph_traits::null_halfedge()) + { + bh = next(bh, g); + if(target(bh, g) == target(ch, g)) + v = ch; + else if(target(bh, g) == target(next(ch, g), g)) + v = next(ch, g); + else if(target(bh, g) == target(next(next(ch, g), g), g)) + v = next(next(ch, g), g); + + if(v == boost::graph_traits::null_halfedge()) + { + bh = next(bh, g); + if(target(bh, g) == target(ch, g)) + v = ch; + else if(target(bh, g) == target(next(ch, g), g)) + v = next(ch, g); + else if(target(bh, g) == target(next(next(ch, g), g), g)) + v = next(next(ch, g), g); + } + } + + return v; + } + void operator()(const Box* b, const Box* c) const { - halfedge_descriptor h = halfedge(b->info(), m_tmesh); - halfedge_descriptor opp_h; + halfedge_descriptor bh = halfedge(b->info(), g); + halfedge_descriptor opp_bh; // check for shared egde - for(unsigned int i=0; i<3; ++i){ - opp_h = opposite(h, m_tmesh); - if(face(opp_h, m_tmesh) == c->info()){ + for(unsigned int i=0; i<3; ++i) + { + opp_bh = opposite(bh, g); + if(face(opp_bh, g) == c->info()) + { // there is an intersection if the four points are coplanar and // the triangles overlap - if(CGAL::coplanar(get(m_vpmap, target(h, m_tmesh)), - get(m_vpmap, target(next(h, m_tmesh), m_tmesh)), - get(m_vpmap, source(h, m_tmesh)), - get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh))) && - CGAL::coplanar_orientation(get(m_vpmap, source(h, m_tmesh)), - get(m_vpmap, target(h, m_tmesh)), - get(m_vpmap, target(next(h, m_tmesh), m_tmesh)), - get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh))) - == CGAL::POSITIVE){ + if(CGAL::coplanar(get(m_vpmap, target(bh, g)), + get(m_vpmap, target(next(bh, g), g)), + get(m_vpmap, source(bh, g)), + get(m_vpmap, target(next(opp_bh, g), g))) && + CGAL::coplanar_orientation(get(m_vpmap, source(bh, g)), + get(m_vpmap, target(bh, g)), + get(m_vpmap, target(next(bh, g), g)), + get(m_vpmap, target(next(opp_bh, g), g))) == CGAL::POSITIVE) + { *m_iterator_wrapper++ = std::make_pair(b->info(), c->info()); return; } else { // there is a shared edge but no intersection return; } } - h = next(h, m_tmesh); + + bh = next(bh, g); } // check for shared vertex --> maybe intersection, maybe not - halfedge_descriptor g = halfedge(c->info(),m_tmesh); - halfedge_descriptor v; - - if(target(h,m_tmesh) == target(g,m_tmesh)) - v = g; - if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh)) - v = next(g,m_tmesh); - if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh)) - v = next(next(g,m_tmesh),m_tmesh); - - if(v == halfedge_descriptor()){ - h = next(h,m_tmesh); - if(target(h,m_tmesh) == target(g,m_tmesh)) - v = g; - if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh)) - v = next(g,m_tmesh); - if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh)) - v = next(next(g,m_tmesh),m_tmesh); - if(v == halfedge_descriptor()){ - h = next(h,m_tmesh); - if(target(h,m_tmesh) == target(g,m_tmesh)) - v = g; - if(target(h,m_tmesh) == target(next(g,m_tmesh),m_tmesh)) - v = next(g,m_tmesh); - if(target(h,m_tmesh) == target(next(next(g,m_tmesh),m_tmesh),m_tmesh)) - v = next(next(g,m_tmesh),m_tmesh); - } - } - - if(v != halfedge_descriptor()){ + halfedge_descriptor ch = halfedge(c->info(), g); + halfedge_descriptor v = common_vertex(bh, ch); + if(v != boost::graph_traits::null_halfedge()) + { // found shared vertex: - CGAL_assertion(target(h,m_tmesh) == target(v,m_tmesh)); + CGAL_assertion(target(bh, g) == target(v, g)); + // geometric check if the opposite segments intersect the triangles - Triangle t1 = triangle_functor( get(m_vpmap,target(h,m_tmesh)), - get(m_vpmap, target(next(h,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh))); - Triangle t2 = triangle_functor( get(m_vpmap, target(v,m_tmesh)), - get(m_vpmap, target(next(v,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(v,m_tmesh),m_tmesh),m_tmesh))); + // Note that there is no shared edge here, the case was handled above. + Triangle t1 = triangle_functor(get(m_vpmap,target(bh, g)), + get(m_vpmap, target(next(bh, g), g)), + get(m_vpmap, target(next(next(bh, g), g), g))); + Triangle t2 = triangle_functor(get(m_vpmap, target(v, g)), + get(m_vpmap, target(next(v, g), g)), + get(m_vpmap, target(next(next(v, g), g), g))); - Segment s1 = segment_functor( get(m_vpmap, target(next(h,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh))); - Segment s2 = segment_functor( get(m_vpmap, target(next(v,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(v,m_tmesh),m_tmesh),m_tmesh))); + Segment s1 = segment_functor(get(m_vpmap, target(next(bh, g), g)), + get(m_vpmap, target(next(next(bh, g), g), g))); + Segment s2 = segment_functor(get(m_vpmap, target(next(v, g), g)), + get(m_vpmap, target(next(next(v, g), g), g))); - if(do_intersect_3_functor(t1,s2)){ + if(do_intersect_3_functor(t1, s2)) *m_iterator_wrapper++ = std::make_pair(b->info(), c->info()); - } else if(do_intersect_3_functor(t2,s1)){ + else if(do_intersect_3_functor(t2, s1)) *m_iterator_wrapper++ = std::make_pair(b->info(), c->info()); - } + return; } // check for geometric intersection - Triangle t1 = triangle_functor( get(m_vpmap, target(h,m_tmesh)), - get(m_vpmap, target(next(h,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(h,m_tmesh),m_tmesh),m_tmesh))); - Triangle t2 = triangle_functor( get(m_vpmap, target(g,m_tmesh)), - get(m_vpmap, target(next(g,m_tmesh),m_tmesh)), - get(m_vpmap, target(next(next(g,m_tmesh),m_tmesh),m_tmesh))); - if(do_intersect_3_functor(t1, t2)){ + Triangle t1 = triangle_functor(get(m_vpmap, target(bh, g)), + get(m_vpmap, target(next(bh, g), g)), + get(m_vpmap, target(next(next(bh, g), g), g))); + Triangle t2 = triangle_functor(get(m_vpmap, target(ch, g)), + get(m_vpmap, target(next(ch, g), g)), + get(m_vpmap, target(next(next(ch, g), g), g))); + + if(do_intersect_3_functor(t1, t2)) *m_iterator_wrapper++ = std::make_pair(b->info(), c->info()); - } - } // end operator () + } }; // end struct Intersect_facets struct Throw_at_output { diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_predicates.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_predicates.h index 6cc49495ad5..d4fdaa2f018 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_predicates.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_predicates.h @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -85,6 +86,81 @@ bool is_degenerate_edge(typename boost::graph_traits::edge_descript return is_degenerate_edge(e, pm, parameters::all_default()); } +/// \ingroup PMP_repairing_grp +/// collects the degenerate edges within a given range of edges. +/// +/// @tparam EdgeRange a model of `Range` with value type `boost::graph_traits::%edge_descriptor` +/// @tparam TriangleMesh a model of `HalfedgeGraph` +/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" +/// +/// @param edges a subset of edges of `tm` +/// @param tm a triangle mesh +/// @param out an output iterator in which the degenerate edges are written +/// @param np optional \ref pmp_namedparameters "Named Parameters" described below +/// +/// \cgalNamedParamsBegin +/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`. +/// The type of this map is model of `ReadWritePropertyMap`. +/// If this parameter is omitted, an internal property map for +/// `CGAL::vertex_point_t` should be available in `TriangleMesh` +/// \cgalParamEnd +/// \cgalParamBegin{geom_traits} a geometric traits class instance. +/// The traits class must provide the nested type `Point_3`, +/// and the nested functor `Equal_3` to check whether two points are identical. +/// \cgalParamEnd +/// \cgalNamedParamsEnd +template +OutputIterator degenerate_edges(const EdgeRange& edges, + const TriangleMesh& tm, + OutputIterator out, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + for(edge_descriptor ed : edges) + if(is_degenerate_edge(ed, tm, np)) + *out++ = ed; + + return out; +} + +template +OutputIterator degenerate_edges(const EdgeRange& edges, + const TriangleMesh& tm, + OutputIterator out, + typename boost::enable_if< + typename boost::has_range_iterator + >::type* = 0) +{ + return degenerate_edges(edges, tm, out, CGAL::parameters::all_default()); +} + +/// \ingroup PMP_repairing_grp +/// calls the function `degenerate_edges()` with the range: `edges(tm)`. +/// +/// See above for the comprehensive description of the parameters. +/// +template +OutputIterator degenerate_edges(const TriangleMesh& tm, + OutputIterator out, + const NamedParameters& np +#ifndef DOXYGEN_RUNNING + , + typename boost::disable_if< + boost::has_range_iterator + >::type* = 0 +#endif + ) +{ + return degenerate_edges(edges(tm), tm, out, np); +} + +template +OutputIterator degenerate_edges(const TriangleMesh& tm, OutputIterator out) +{ + return degenerate_edges(edges(tm), tm, out, CGAL::parameters::all_default()); +} + /// \ingroup PMP_repairing_grp /// checks whether a triangle face is degenerate. /// A triangle face is considered degenerate if the geometric positions of its vertices are collinear. @@ -142,6 +218,83 @@ bool is_degenerate_triangle_face(typename boost::graph_traits::fac return CGAL::Polygon_mesh_processing::is_degenerate_triangle_face(f, tm, parameters::all_default()); } +/// \ingroup PMP_repairing_grp +/// collects the degenerate faces within a given range of faces. +/// +/// @tparam FaceRange a model of `Range` with value type `boost::graph_traits::%face_descriptor` +/// @tparam TriangleMesh a model of `FaceGraph` +/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters" +/// +/// @param faces a subset of faces of `tm` +/// @param tm a triangle mesh +/// @param out an output iterator in which the degenerate faces are put +/// @param np optional \ref pmp_namedparameters "Named Parameters" described below +/// +/// \cgalNamedParamsBegin +/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`. +/// The type of this map is model of `ReadWritePropertyMap`. +/// If this parameter is omitted, an internal property map for +/// `CGAL::vertex_point_t` should be available in `TriangleMesh` +/// \cgalParamEnd +/// \cgalParamBegin{geom_traits} a geometric traits class instance. +/// The traits class must provide the nested functor `Collinear_3` +/// to check whether three points are collinear. +/// \cgalParamEnd +/// \cgalNamedParamsEnd +/// +template +OutputIterator degenerate_faces(const FaceRange& faces, + const TriangleMesh& tm, + OutputIterator out, + const NamedParameters& np) +{ + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + for(face_descriptor fd : faces) + { + if(is_degenerate_triangle_face(fd, tm, np)) + *out++ = fd; + } + return out; +} + +template +OutputIterator degenerate_faces(const FaceRange& faces, + const TriangleMesh& tm, + OutputIterator out, + typename boost::enable_if< + boost::has_range_iterator + >::type* = 0) +{ + return degenerate_faces(faces, tm, out, CGAL::parameters::all_default()); +} + +/// \ingroup PMP_repairing_grp +/// calls the function `degenerate_faces()` with the range: `faces(tm)`. +/// +/// See above for the comprehensive description of the parameters. +/// +template +OutputIterator degenerate_faces(const TriangleMesh& tm, + OutputIterator out, + const NamedParameters& np +#ifndef DOXYGEN_RUNNING + , + typename boost::disable_if< + boost::has_range_iterator + >::type* = 0 +#endif + ) +{ + return degenerate_faces(faces(tm), tm, out, np); +} + +template +OutputIterator degenerate_faces(const TriangleMesh& tm, OutputIterator out) +{ + return degenerate_faces(faces(tm), tm, out, CGAL::parameters::all_default()); +} + /// \ingroup PMP_repairing_grp /// checks whether a triangle face is needle. /// A triangle is said to be a needle if its longest edge is much longer than its shortest edge. diff --git a/Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h b/Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h index 82f712fdc16..f3f01579d2c 100644 --- a/Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h +++ b/Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h @@ -32,6 +32,8 @@ #include #include #include +#include +#include #include #include #include @@ -45,7 +47,7 @@ #include #include #include -#include +#include // the named parameter header being not documented the doc is put here for now #ifdef DOXYGEN_RUNNING diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 2bf92f35e85..c047f553df6 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -98,6 +98,7 @@ endif() create_single_source_cgal_program("test_pmp_manifoldness.cpp") create_single_source_cgal_program("test_mesh_smoothing.cpp") create_single_source_cgal_program("test_remove_caps_needles.cpp") + create_single_source_cgal_program("test_pmp_remove_self_intersections.cpp") if( TBB_FOUND ) CGAL_target_use_TBB(test_pmp_distance) @@ -114,6 +115,8 @@ find_package(Ceres QUIET) if(TARGET ceres) target_compile_definitions( test_mesh_smoothing PRIVATE CGAL_PMP_USE_CERES_SOLVER ) target_link_libraries( test_mesh_smoothing PRIVATE ceres ) + target_compile_definitions( test_pmp_remove_self_intersections PRIVATE CGAL_PMP_USE_CERES_SOLVER ) + target_link_libraries( test_pmp_remove_self_intersections PRIVATE ceres ) endif(TARGET ceres) if(BUILD_TESTING) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_remove_self_intersections.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_remove_self_intersections.cpp new file mode 100644 index 00000000000..e84c9120009 --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_remove_self_intersections.cpp @@ -0,0 +1,126 @@ +#include +#include + +#include +#include + +#include +#include +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +typedef CGAL::Surface_mesh Surface_mesh; + +typedef boost::graph_traits::edge_descriptor edge_descriptor; +typedef boost::graph_traits::face_descriptor face_descriptor; + +void fix_self_intersections(const char* fname) +{ + std::cout << std::endl << "---------------" << std::endl; + std::cout << "Test " << fname << std::endl; + + std::ifstream input(fname); + Surface_mesh mesh; + if (!input || !(input >> mesh) || mesh.is_empty()) { + std::cerr << "Error: " << fname << " is not a valid off file.\n"; + std::exit(1); + } + + Surface_mesh mesh_cpy = mesh; + + PMP::remove_self_intersections(mesh); + assert( CGAL::is_valid_polygon_mesh(mesh) ); + CGAL_warning( !PMP::does_self_intersect(mesh) ); + + std::ofstream out("post_repair.off"); + out.precision(17); + out << mesh; + out.close(); + + // just to check compilation + PMP::remove_self_intersections(mesh, CGAL::parameters::number_of_iterations(10)); +} + +void fix_local_self_intersections(const char* mesh_filename, const char* mesh_selection_filename) +{ + std::cout << std::endl << "---------------" << std::endl; + std::cout << "Test " << mesh_filename << " with selection " << mesh_selection_filename << std::endl; + + std::ifstream input(mesh_filename); + Surface_mesh mesh; + if (!input || !(input >> mesh) || mesh.is_empty()) { + std::cerr << "Error: " << mesh_filename << " is not a valid off file.\n"; + std::exit(1); + } + + std::ifstream selection_input(mesh_selection_filename); + std::list selected_faces; + std::string line; + // skip the first line (faces are on the second line) + if(!selection_input || !std::getline(selection_input, line) || !std::getline(selection_input, line)) + { + std::cerr << "Error: could not read selection: " << mesh_selection_filename << std::endl; + std::exit(1); + } + + std::istringstream face_line(line); + std::size_t face_id; + while(face_line >> face_id) + selected_faces.push_back(*(faces(mesh).begin() + face_id)); + std::cout << selected_faces.size() << " faces selected" << std::endl; + + PMP::remove_self_intersections(selected_faces, mesh, CGAL::parameters::verbosity_level(true)); + assert( CGAL::is_valid_polygon_mesh(mesh) ); + CGAL_warning( !PMP::does_self_intersect(selected_faces, mesh) ); + + std::ofstream out("post_repair.off"); + out.precision(17); + out << mesh; + out.close(); + + std::unordered_map face_index_map; + int counter = 0; + BOOST_FOREACH(face_descriptor fd, faces(mesh)) + face_index_map[fd] = counter++; + + std::ofstream out_final_selection("post_repair.selection.txt"); + out_final_selection << std::endl; // first line are vertex indices + + BOOST_FOREACH(const face_descriptor fd, selected_faces) + out_final_selection << face_index_map[fd] << " "; + out_final_selection << std::endl << std::endl; + + PMP::remove_self_intersections(selected_faces, mesh); + assert( CGAL::is_valid_polygon_mesh(mesh) ); + CGAL_warning( !PMP::does_self_intersect(selected_faces, mesh) ); +} + +int main() +{ +#if 0 + fix_self_intersections("data_repair/brain.off"); + fix_self_intersections("data_repair/flute.off"); + fix_self_intersections("data_repair/dinosaur.off"); + fix_self_intersections("data_repair/hemispheres.off"); + + // selection is adjacent to a self-intersection but does not contain any intersection + fix_local_self_intersections("data_repair/brain.off", "data_repair/brain-complete.selection.txt"); + + // selection covers nicely a self-intersection + fix_local_self_intersections("data_repair/brain.off", "data_repair/brain-adjacent.selection.txt"); + + // selection contains part of a self intersection + fix_local_self_intersections("data_repair/brain.off", "data_repair/brain-partial.selection.txt"); + + // selection contains disjoint parts of a self intersection + fix_local_self_intersections("data_repair/brain.off", "data_repair/brain-disjoint.selection.txt"); +#endif + + // Remove only self-intersections within a single hemisphere + fix_local_self_intersections("data_repair/hemispheres.off", "data_repair/hemispheres-half.selection.txt"); + + return 0; +}