Merge branch 'PMP-Make_remove_self_intersections_local-GF-old' into PMP-Make_remove_self_intersections_local-GF

This commit is contained in:
Mael Rouxel-Labbé 2020-02-03 13:05:52 +01:00
commit 32aa902bf5
13 changed files with 4527 additions and 3024 deletions

View File

@ -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)

View File

@ -0,0 +1,83 @@
#define CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/remove_self_intersections.h>
#include <CGAL/Polygon_mesh_processing/remove_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/IO/STL_reader.h>
#include <fstream>
#include <sstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef typename boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Mesh>::face_descriptor face_descriptor;
namespace PMP = CGAL::Polygon_mesh_processing;
template <typename K, typename Mesh>
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<Point> points;
std::vector<std::vector<int> > 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<K>(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;
}

View File

@ -566,10 +566,8 @@ double approximate_Hausdorff_distance(
VertexPointMap vpm_2)
{
std::vector<typename Kernel::Point_3> 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<Concurrency_tag, Kernel>(sample_points, tm2, vpm_2);
}

View File

@ -21,7 +21,7 @@
#endif
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
@ -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<edge_descriptor, bool> 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<edge_descriptor, bool> other_hd_already_exists = edge(v2, v3, mesh_);
if(other_hd_already_exists.second)
return false;
return true;
return (alpha + beta > CGAL_PI);
}
template <typename Marked_edges_map, typename EdgeRange>

View File

@ -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 <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/algorithm.h>
#include <CGAL/assertions.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/iterator.h>
#include <CGAL/property_map.h>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
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 <typename PolygonMesh>
bool is_non_manifold_vertex(typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
const PolygonMesh& pm)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef CGAL::dynamic_halfedge_property_t<bool> Halfedge_property_tag;
typedef typename boost::property_map<PolygonMesh, Halfedge_property_tag>::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 <typename G>
struct Vertex_collector
{
typedef typename boost::graph_traits<G>::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<vertex_descriptor>& verts = collections[v1];
if(verts.empty())
verts.push_back(v1);
verts.push_back(v2);
}
template<typename OutputIterator>
void dump(OutputIterator out)
{
typedef std::pair<const vertex_descriptor, std::vector<vertex_descriptor> > Pair_type;
for(const Pair_type& p : collections)
*out++ = p.second;
}
void dump(Emptyset_iterator) { }
std::map<vertex_descriptor, std::vector<vertex_descriptor> > collections;
};
template <typename PolygonMesh, typename VPM, typename ConstraintMap>
typename boost::graph_traits<PolygonMesh>::vertex_descriptor
create_new_vertex_for_sector(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor sector_begin_h,
typename boost::graph_traits<PolygonMesh>::halfedge_descriptor sector_last_h,
PolygonMesh& pm,
const VPM& vpm,
const ConstraintMap& cmap)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::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 <typename PolygonMesh, typename NamedParameters>
std::size_t make_umbrella_manifold(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
PolygonMesh& pm,
internal::Vertex_collector<PolygonMesh>& dmap,
const NamedParameters& np)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
using parameters::get_parameter;
using parameters::choose_parameter;
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::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<internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool> // default (no constraint pmap)
>::type VerticesMap;
VerticesMap cmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(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 <i>umbrella</i> 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<PolygonMesh>::%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 <typename PolygonMesh, typename OutputIterator>
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<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef CGAL::dynamic_vertex_property_t<bool> Vertex_bool_tag;
typedef typename boost::property_map<PolygonMesh, Vertex_bool_tag>::const_type Known_manifold_vertex_map;
typedef CGAL::dynamic_vertex_property_t<halfedge_descriptor> Vertex_halfedge_tag;
typedef typename boost::property_map<PolygonMesh, Vertex_halfedge_tag>::const_type Visited_vertex_map;
typedef CGAL::dynamic_halfedge_property_t<bool> Halfedge_property_tag;
typedef typename boost::property_map<PolygonMesh, Halfedge_property_tag>::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<PolygonMesh>::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<vertex_descriptor>`. 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 <typename PolygonMesh, typename NamedParameters>
std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm,
const NamedParameters& np)
{
using parameters::get_parameter;
using parameters::choose_parameter;
typedef boost::graph_traits<PolygonMesh> 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<halfedge_descriptor> non_manifold_cones;
non_manifold_vertices(pm, std::back_inserter(non_manifold_cones));
internal::Vertex_collector<PolygonMesh> 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 <class PolygonMesh>
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

View File

@ -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<TM, boost::vertex_point_t>::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<TM>::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<TM>::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<TM>::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<TM>::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 {

View File

@ -23,6 +23,7 @@
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/helpers.h>
#include <boost/range/has_range_iterator.hpp>
#include <boost/graph/graph_traits.hpp>
#include <limits>
@ -85,6 +86,81 @@ bool is_degenerate_edge(typename boost::graph_traits<PolygonMesh>::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<TriangleMesh>::%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 <typename EdgeRange, typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_edges(const EdgeRange& edges,
const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np)
{
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
for(edge_descriptor ed : edges)
if(is_degenerate_edge(ed, tm, np))
*out++ = ed;
return out;
}
template <typename EdgeRange, typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_edges(const EdgeRange& edges,
const TriangleMesh& tm,
OutputIterator out,
typename boost::enable_if<
typename boost::has_range_iterator<EdgeRange>
>::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 <typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_edges(const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np
#ifndef DOXYGEN_RUNNING
,
typename boost::disable_if<
boost::has_range_iterator<TriangleMesh>
>::type* = 0
#endif
)
{
return degenerate_edges(edges(tm), tm, out, np);
}
template <typename TriangleMesh, typename OutputIterator>
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<TriangleMesh>::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<TriangleMesh>::%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 <typename FaceRange, typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_faces(const FaceRange& faces,
const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np)
{
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
for(face_descriptor fd : faces)
{
if(is_degenerate_triangle_face(fd, tm, np))
*out++ = fd;
}
return out;
}
template <typename FaceRange, typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_faces(const FaceRange& faces,
const TriangleMesh& tm,
OutputIterator out,
typename boost::enable_if<
boost::has_range_iterator<FaceRange>
>::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 <typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_faces(const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np
#ifndef DOXYGEN_RUNNING
,
typename boost::disable_if<
boost::has_range_iterator<TriangleMesh>
>::type* = 0
#endif
)
{
return degenerate_faces(faces(tm), tm, out, np);
}
template <typename TriangleMesh, typename OutputIterator>
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 <i>needle</i> if its longest edge is much longer than its shortest edge.

View File

@ -32,6 +32,8 @@
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/remove_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/remove_self_intersections.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/detect_features.h>
@ -45,7 +47,7 @@
#include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include <CGAL/Polygon_mesh_processing/internal/remove_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/manifoldness.h>
// the named parameter header being not documented the doc is put here for now
#ifdef DOXYGEN_RUNNING

View File

@ -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)

View File

@ -0,0 +1,126 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <iostream>
#include <fstream>
#include <unordered_map>
namespace PMP = CGAL::Polygon_mesh_processing;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<Surface_mesh>::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<face_descriptor> 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_descriptor, int> 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;
}