songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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-07 17:09:59 +01:00
parent 9892729a20 119f5b913b
commit cb3e5cd83c
26 changed files with 5379 additions and 3786 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Kernel_23/doc/Kernel_23/CGAL/Vector_3.h
View File

@ -68,12 +68,12 @@ introduces a vector `v` initialized to `(x, y, z)`.
Vector_3(int x, int y, int z); Vector_3(int x, int y, int z);
/*! /*!
introduces a vector `v` initialized to `(x, y, z). introduces a vector `v` initialized to `(x, y, z)`.
*/ */
Vector_3(double x, double y, double z); Vector_3(double x, double y, double z);
/*! /*!
introduces a vector `v` initialized to `(hx/hw, hy/hw, hz/hw). introduces a vector `v` initialized to `(hx/hw, hy/hw, hz/hw)`.
*/ */
Vector_3(const Kernel::RT &hx, const Kernel::RT &hy, const Kernel::RT &hz, const Kernel::RT &hw = RT(1)); Vector_3(const Kernel::RT &hx, const Kernel::RT &hy, const Kernel::RT &hz, const Kernel::RT &hw = RT(1));

4
Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt
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( "repair_polygon_soup_example.cpp" )
create_single_source_cgal_program( "mesh_smoothing_example.cpp") create_single_source_cgal_program( "mesh_smoothing_example.cpp")
create_single_source_cgal_program( "locate_example.cpp") create_single_source_cgal_program( "locate_example.cpp")
create_single_source_cgal_program( "mesh_self_intersection_removal_example.cpp")
if(OpenMesh_FOUND) if(OpenMesh_FOUND)
create_single_source_cgal_program( "compute_normals_example_OM.cpp" ) create_single_source_cgal_program( "compute_normals_example_OM.cpp" )
@ -119,5 +120,8 @@ find_package(Ceres QUIET)
if(TARGET ceres) if(TARGET ceres)
target_compile_definitions( mesh_smoothing_example PRIVATE CGAL_PMP_USE_CERES_SOLVER ) target_compile_definitions( mesh_smoothing_example PRIVATE CGAL_PMP_USE_CERES_SOLVER )
target_link_libraries( mesh_smoothing_example PRIVATE ceres ) 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) endif(TARGET ceres)

83
Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_self_intersection_removal_example.cpp Normal file
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::experimental::remove_self_intersections(mesh);
std::ofstream("out.off") << std::setprecision(17) << mesh;
std::cout << "Success? " << !PMP::does_self_intersect(mesh) << std::endl;
return 0;
}

24
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h
View File

@ -98,6 +98,12 @@ void sum_normals(const PM& pmesh,
const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh)); const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh));
const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits); const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits);
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Normal of " << f << " pts: " << pv << " ; " << pvn << " ; " << pvnn << std::endl;
std::cout << " --> " << n << std::endl;
#endif
sum = traits.construct_sum_of_vectors_3_object()(sum, n); sum = traits.construct_sum_of_vectors_3_object()(sum, n);
he = next(he, pmesh); he = next(he, pmesh);
@ -206,7 +212,7 @@ void compute_face_normals(const PolygonMesh& pmesh,
{ {
typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np); typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np);
put(face_normals, f, vec); put(face_normals, f, vec);
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "normal at face " << f << " is " << get(face_normals, f) << std::endl; std::cout << "normal at face " << f << " is " << get(face_normals, f) << std::endl;
#endif #endif
} }
@ -522,6 +528,10 @@ compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<Po
typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object(); typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object();
typename GT::Compute_squared_length_3 csl_3 = traits.compute_squared_length_3_object(); typename GT::Compute_squared_length_3 csl_3 = traits.compute_squared_length_3_object();
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Compute normal as weighted sums; type: " << vn_type << std::endl;
#endif
Vector_3 normal = cv_3(CGAL::NULL_VECTOR); Vector_3 normal = cv_3(CGAL::NULL_VECTOR);
halfedge_descriptor h = halfedge(v, pmesh); halfedge_descriptor h = halfedge(v, pmesh);
@ -548,7 +558,13 @@ compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<Po
const FT den = CGAL::approximate_sqrt(csl_3(v1) * csl_3(v2)); const FT den = CGAL::approximate_sqrt(csl_3(v1) * csl_3(v2));
if(den == FT(0)) if(den == FT(0))
{
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Null denominator, switching to no weights" << std::endl;
#endif
return compute_vertex_normal_as_sum_of_weighted_normals(v, NO_WEIGHT, face_normals, vpmap, pmesh, traits); return compute_vertex_normal_as_sum_of_weighted_normals(v, NO_WEIGHT, face_normals, vpmap, pmesh, traits);
}
n = traits.construct_scaled_vector_3_object()(n, FT(1) / den); n = traits.construct_scaled_vector_3_object()(n, FT(1) / den);
normal = traits.construct_sum_of_vectors_3_object()(normal, n); normal = traits.construct_sum_of_vectors_3_object()(normal, n);
@ -632,9 +648,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal)); const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal));
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << std::endl << std::endl; std::cout << "<----- compute vertex normal at " << get(vpmap, v)
std::cout << "----------------------------------------------------------------------" << std::endl;
std::cout << "compute vertex at " << get(vpmap, v)
<< ", must compute face normals? " << must_compute_face_normals << std::endl; << ", must compute face normals? " << must_compute_face_normals << std::endl;
#endif #endif
@ -654,7 +668,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
} }
} }
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Incident face normals:" << std::endl; std::cout << "Incident face normals:" << std::endl;
for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh)) for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh))
{ {

6
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h
View File

@ -566,10 +566,8 @@ double approximate_Hausdorff_distance(
VertexPointMap vpm_2) VertexPointMap vpm_2)
{ {
std::vector<typename Kernel::Point_3> sample_points; std::vector<typename Kernel::Point_3> sample_points;
sample_triangle_mesh( sample_triangle_mesh(tm1,std::back_inserter(sample_points),np);
tm1,
std::back_inserter(sample_points),
np);
return approximate_Hausdorff_distance<Concurrency_tag, Kernel>(sample_points, tm2, vpm_2); return approximate_Hausdorff_distance<Concurrency_tag, Kernel>(sample_points, tm2, vpm_2);
} }

4
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h
View File

@ -170,7 +170,7 @@ public:
// calls compute once to factorize with the preconditioner // calls compute once to factorize with the preconditioner
if(!solver.factor(A, D)) if(!solver.factor(A, D))
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Could not factorize linear system with preconditioner." << std::endl; std::cerr << "Could not factorize linear system with preconditioner." << std::endl;
#endif #endif
return false; return false;
@ -180,7 +180,7 @@ public:
!solver.linear_solver(by, Xy) || !solver.linear_solver(by, Xy) ||
!solver.linear_solver(bz, Xz)) !solver.linear_solver(bz, Xz))
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Could not solve linear system." << std::endl; std::cerr << "Could not solve linear system." << std::endl;
#endif #endif
return false; return false;

192
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h
View File

@ -21,7 +21,7 @@
#endif #endif
#include <CGAL/Polygon_mesh_processing/compute_normal.h> #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_tree.h>
#include <CGAL/AABB_traits.h> #include <CGAL/AABB_traits.h>
@ -50,9 +50,11 @@ namespace Polygon_mesh_processing {
namespace internal { namespace internal {
template <typename V, typename GT> template <typename V, typename GT>
double get_radian_angle(const V& v1, const V& v2, const GT& gt) typename GT::FT get_radian_angle(const V& v1, const V& v2, const GT& gt)
{ {
return gt.compute_approximate_angle_3_object()(v1, v2) * CGAL_PI / 180.; typedef typename GT::FT FT;
return gt.compute_approximate_angle_3_object()(v1, v2) * CGAL_PI / FT(180);
} }
// super naive for now. Not sure it even makes sense to do something like that for surfaces // super naive for now. Not sure it even makes sense to do something like that for surfaces
@ -68,6 +70,7 @@ class Delaunay_edge_flipper
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor; typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref; typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector; typedef typename GeomTraits::Vector_3 Vector;
public: public:
@ -86,20 +89,13 @@ public:
const halfedge_descriptor h = halfedge(e, mesh_); const halfedge_descriptor h = halfedge(e, mesh_);
const halfedge_descriptor opp_h = opposite(h, mesh_); const halfedge_descriptor opp_h = opposite(h, mesh_);
vertex_descriptor v0 = source(h, mesh_); const vertex_descriptor v0 = source(h, mesh_);
vertex_descriptor v1 = target(h, mesh_); const vertex_descriptor v1 = target(h, mesh_);
vertex_descriptor v2 = target(next(h, mesh_), mesh_); const vertex_descriptor v2 = target(next(h, mesh_), mesh_);
vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_); const vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_);
const Point_ref p0 = get(vpmap_, v0);
const Point_ref p1 = get(vpmap_, v1);
const Point_ref p2 = get(vpmap_, v2);
const Point_ref p3 = get(vpmap_, v3);
double alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_); std::set<vertex_descriptor> unique_vs { v0, v1, v2, v3 };
double beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_); if(unique_vs.size() != 4)
// not local Delaunay if the sum of the angles is greater than pi
if(alpha + beta <= CGAL_PI)
return false; return false;
// Don't want to flip if the other diagonal already exists // Don't want to flip if the other diagonal already exists
@ -108,7 +104,16 @@ public:
if(other_hd_already_exists.second) if(other_hd_already_exists.second)
return false; return false;
return true; // 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);
const Point_ref p3 = get(vpmap_, v3);
FT alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_);
FT beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_);
return (alpha + beta > CGAL_PI);
} }
template <typename Marked_edges_map, typename EdgeRange> template <typename Marked_edges_map, typename EdgeRange>
@ -127,6 +132,10 @@ public:
template <typename FaceRange> template <typename FaceRange>
void operator()(const FaceRange& face_range) void operator()(const FaceRange& face_range)
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Flipping edges" << std::endl;
#endif
// edges to consider // edges to consider
std::vector<edge_descriptor> edge_range; std::vector<edge_descriptor> edge_range;
edge_range.reserve(3 * face_range.size()); edge_range.reserve(3 * face_range.size());
@ -166,6 +175,10 @@ public:
++flipped_n; ++flipped_n;
halfedge_descriptor h = halfedge(e, mesh_); halfedge_descriptor h = halfedge(e, mesh_);
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Flipping " << edge(h, mesh_) << std::endl;
#endif
Euler::flip_edge(h, mesh_); Euler::flip_edge(h, mesh_);
add_to_stack_if_unmarked(edge(next(h, mesh_), mesh_), marks, edge_range); add_to_stack_if_unmarked(edge(next(h, mesh_), mesh_), marks, edge_range);
@ -174,6 +187,10 @@ public:
add_to_stack_if_unmarked(edge(prev(opposite(h, mesh_), mesh_), mesh_), marks, edge_range); add_to_stack_if_unmarked(edge(prev(opposite(h, mesh_), mesh_), mesh_), marks, edge_range);
} }
} }
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << flipped_n << " flips" << std::endl;
#endif
} }
private: private:
@ -190,6 +207,7 @@ class Angle_smoother
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref; typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector; typedef typename GeomTraits::Vector_3 Vector;
typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair; typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair;
@ -231,7 +249,7 @@ public:
Vector operator()(const vertex_descriptor v) const Vector operator()(const vertex_descriptor v) const
{ {
Vector move = CGAL::NULL_VECTOR; Vector move = CGAL::NULL_VECTOR;
double weights_sum = 0.; FT weights_sum = FT(0);
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_)) for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{ {
@ -249,26 +267,30 @@ public:
Vector right_v(pt, right_pt); Vector right_v(pt, right_pt);
// rotate // rotate
double angle = get_radian_angle(right_v, left_v, traits_);
CGAL_warning(angle != 0.); // no degenerate faces is a precondition
if(angle == 0.)
continue;
Vector bisector = rotate_edge(main_he, incident_pair); Vector bisector = rotate_edge(main_he, incident_pair);
double scaling_factor = CGAL::approximate_sqrt( FT scaling_factor = CGAL::approximate_sqrt(
traits_.compute_squared_distance_3_object()(get(vpmap_, source(main_he, mesh_)), traits_.compute_squared_distance_3_object()(get(vpmap_, source(main_he, mesh_)),
get(vpmap_, target(main_he, mesh_)))); get(vpmap_, target(main_he, mesh_))));
bisector = traits_.construct_scaled_vector_3_object()(bisector, scaling_factor); bisector = traits_.construct_scaled_vector_3_object()(bisector, scaling_factor);
Vector ps_psi(ps, traits_.construct_translated_point_3_object()(pt, bisector)); Vector ps_psi(ps, traits_.construct_translated_point_3_object()(pt, bisector));
FT angle = get_radian_angle(right_v, left_v, traits_);
if(angle == FT(0))
{
// no degenerate faces is a precondition, angle can be 0 but it should be a numerical error
CGAL_warning(!is_degenerate_triangle_face(face(main_he, mesh_), mesh_));
return ps_psi; // since a small angle gives more weight, a null angle give priority (?)
}
// small angles carry more weight // small angles carry more weight
double weight = 1. / (angle*angle); FT weight = 1. / CGAL::square(angle);
weights_sum += weight; weights_sum += weight;
move += weight * ps_psi; move += weight * ps_psi;
} }
if(weights_sum != 0.) if(weights_sum != FT(0))
move /= weights_sum; move /= weights_sum;
return move; return move;
@ -288,6 +310,7 @@ class Area_smoother
typedef typename boost::property_traits<VertexPointMap>::value_type Point; typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref; typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector; typedef typename GeomTraits::Vector_3 Vector;
public: public:
@ -298,27 +321,27 @@ public:
{ } { }
private: private:
double element_area(const vertex_descriptor v1, FT element_area(const vertex_descriptor v1,
const vertex_descriptor v2, const vertex_descriptor v2,
const vertex_descriptor v3) const const vertex_descriptor v3) const
{ {
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(get(vpmap_, v1), return CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(get(vpmap_, v1),
get(vpmap_, v2), get(vpmap_, v2),
get(vpmap_, v3)))); get(vpmap_, v3)));
} }
double element_area(const Point& P, FT element_area(const Point& P,
const vertex_descriptor v2, const vertex_descriptor v2,
const vertex_descriptor v3) const const vertex_descriptor v3) const
{ {
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(P, return CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(P,
get(vpmap_, v2), get(vpmap_, v2),
get(vpmap_, v3)))); get(vpmap_, v3)));
} }
double compute_average_area_around(const vertex_descriptor v) const FT compute_average_area_around(const vertex_descriptor v) const
{ {
double sum_areas = 0.; FT sum_areas = 0;
unsigned int number_of_edges = 0; unsigned int number_of_edges = 0;
for(halfedge_descriptor h : halfedges_around_source(v, mesh_)) for(halfedge_descriptor h : halfedges_around_source(v, mesh_))
@ -327,7 +350,7 @@ private:
vertex_descriptor vi = source(next(h, mesh_), mesh_); vertex_descriptor vi = source(next(h, mesh_), mesh_);
vertex_descriptor vj = target(next(h, mesh_), mesh_); vertex_descriptor vj = target(next(h, mesh_), mesh_);
double S = element_area(v, vi, vj); FT S = element_area(v, vi, vj);
sum_areas += S; sum_areas += S;
++number_of_edges; ++number_of_edges;
} }
@ -337,7 +360,7 @@ private:
struct Face_energy struct Face_energy
{ {
Face_energy(const Point& pi, const Point& pj, const double s_av) Face_energy(const Point& pi, const Point& pj, const FT s_av)
: :
qx(pi.x()), qy(pi.y()), qz(pi.z()), qx(pi.x()), qy(pi.y()), qz(pi.z()),
rx(pj.x()), ry(pj.y()), rz(pj.z()), rx(pj.x()), ry(pj.y()), rz(pj.z()),
@ -346,7 +369,7 @@ private:
// next two functions are just for convencience, the only thing ceres cares about is the operator() // next two functions are just for convencience, the only thing ceres cares about is the operator()
template <typename T> template <typename T>
double area(const T x, const T y, const T z) const FT area(const T x, const T y, const T z) const
{ {
return CGAL::approximate_sqrt(CGAL::squared_area(Point(x, y, z), return CGAL::approximate_sqrt(CGAL::squared_area(Point(x, y, z),
Point(qx, qy, qz), Point(qx, qy, qz),
@ -354,7 +377,7 @@ private:
} }
template <typename T> template <typename T>
double evaluate(const T x, const T y, const T z) const { return area(x, y, z) - s_av; } FT evaluate(const T x, const T y, const T z) const { return area(x, y, z) - s_av; }
template <typename T> template <typename T>
bool operator()(const T* const x, const T* const y, const T* const z, bool operator()(const T* const x, const T* const y, const T* const z,
@ -390,9 +413,9 @@ private:
} }
private: private:
const double qx, qy, qz; const FT qx, qy, qz;
const double rx, ry, rz; const FT rx, ry, rz;
const double s_av; const FT s_av;
}; };
public: public:
@ -401,12 +424,12 @@ public:
#ifdef CGAL_PMP_USE_CERES_SOLVER #ifdef CGAL_PMP_USE_CERES_SOLVER
const Point_ref vp = get(vpmap_, v); const Point_ref vp = get(vpmap_, v);
const double S_av = compute_average_area_around(v); const FT S_av = compute_average_area_around(v);
const double initial_x = vp.x(); const FT initial_x = vp.x();
const double initial_y = vp.y(); const FT initial_y = vp.y();
const double initial_z = vp.z(); const FT initial_z = vp.z();
double x = initial_x, y = initial_y, z = initial_z; FT x = initial_x, y = initial_y, z = initial_z;
ceres::Problem problem; ceres::Problem problem;
@ -523,7 +546,8 @@ public:
Optimizer compute_move(mesh_, vpmap_, traits_); Optimizer compute_move(mesh_, vpmap_, traits_);
#ifdef CGAL_PMP_SMOOTHING_DEBUG #ifdef CGAL_PMP_SMOOTHING_DEBUG
double total_displacement = 0; FT total_displacement = 0;
std::cout << "apply_moves_in_single_batch: " << apply_moves_in_single_batch << std::endl;
#endif #endif
std::size_t moved_points = 0; std::size_t moved_points = 0;
@ -532,6 +556,10 @@ public:
if(is_border(v, mesh_) || is_constrained(v)) if(is_border(v, mesh_) || is_constrained(v))
continue; continue;
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Considering " << v << " pos: " << get(vpmap_, v) << std::endl;
#endif
// compute normal to v // compute normal to v
Vector vn = compute_vertex_normal(v, mesh_, CGAL::parameters::vertex_point_map(vpmap_) Vector vn = compute_vertex_normal(v, mesh_, CGAL::parameters::vertex_point_map(vpmap_)
.geom_traits(traits_)); .geom_traits(traits_));
@ -542,17 +570,17 @@ public:
// Gram Schmidt so that the new location is on the tangent plane of v (i.e. do mv -= (mv*n)*n) // Gram Schmidt so that the new location is on the tangent plane of v (i.e. do mv -= (mv*n)*n)
const FT sp = traits_.compute_scalar_product_3_object()(vn, move); const FT sp = traits_.compute_scalar_product_3_object()(vn, move);
move = traits_.construct_sum_of_vectors_3_object()(move, move = traits_.construct_sum_of_vectors_3_object()(
traits_.construct_scaled_vector_3_object()(vn, - sp)); move, traits_.construct_scaled_vector_3_object()(vn, - sp));
Point new_pos = pos + move; const Point new_pos = pos + move;
if(move != CGAL::NULL_VECTOR &&
if((!use_sanity_checks || !does_move_create_bad_faces(v, new_pos)) && !does_move_create_degenerate_faces(v, new_pos) &&
(!use_sanity_checks || !does_move_create_bad_faces(v, new_pos)) &&
(!enforce_no_min_angle_regression || does_improve_min_angle_in_star(v, new_pos))) (!enforce_no_min_angle_regression || does_improve_min_angle_in_star(v, new_pos)))
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG #ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "moving " << get(vpmap_, v) << " to " << new_pos << std::endl; std::cout << "moving " << get(vpmap_, v) << " to " << new_pos << std::endl;
total_displacement += CGAL::approximate_sqrt(traits_.compute_squared_length_3_object()(move));
#endif #endif
if(apply_moves_in_single_batch) if(apply_moves_in_single_batch)
@ -560,11 +588,15 @@ public:
else else
put(vpmap_, v, new_pos); put(vpmap_, v, new_pos);
#ifdef CGAL_PMP_SMOOTHING_DEBUG
total_displacement += CGAL::approximate_sqrt(traits_.compute_squared_length_3_object()(move));
#endif
++moved_points; ++moved_points;
} }
else // some sanity check failed else // some sanity check failed
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG #ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "move is rejected!" << std::endl; std::cout << "move is rejected!" << std::endl;
#endif #endif
if(apply_moves_in_single_batch) if(apply_moves_in_single_batch)
@ -607,6 +639,10 @@ public:
Point_ref p_query = get(vpmap_, v); Point_ref p_query = get(vpmap_, v);
const Point projected = tree.closest_point(p_query); const Point projected = tree.closest_point(p_query);
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << p_query << " to " << projected << std::endl;
#endif
put(vpmap_, v, projected); put(vpmap_, v, projected);
} }
} }
@ -617,11 +653,10 @@ private:
return get(vcmap_, v); return get(vcmap_, v);
} }
// check for degenerate or inversed faces // Null faces are bad because they make normal computation difficult
bool does_move_create_bad_faces(const vertex_descriptor v, bool does_move_create_degenerate_faces(const vertex_descriptor v,
const Point& new_pos) const const Point& new_pos) const
{ {
// check for null faces and face inversions
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_)) for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{ {
const halfedge_descriptor prev_he = prev(main_he, mesh_); const halfedge_descriptor prev_he = prev(main_he, mesh_);
@ -630,6 +665,23 @@ private:
if(traits_.collinear_3_object()(lpt, rpt, new_pos)) if(traits_.collinear_3_object()(lpt, rpt, new_pos))
return true; return true;
}
return false;
}
// check for degenerate or inversed faces
bool does_move_create_bad_faces(const vertex_descriptor v,
const Point& new_pos) const
{
// check for face inversions
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const halfedge_descriptor prev_he = prev(main_he, mesh_);
const Point_ref lpt = get(vpmap_, target(main_he, mesh_));
const Point_ref rpt = get(vpmap_, source(prev_he, mesh_));
CGAL_assertion(!traits_.collinear_3_object()(lpt, rpt, new_pos)); // checked above
const Point_ref old_pos = get(vpmap_, v); const Point_ref old_pos = get(vpmap_, v);
Vector ov_1 = traits_.construct_vector_3_object()(old_pos, lpt); Vector ov_1 = traits_.construct_vector_3_object()(old_pos, lpt);
@ -641,7 +693,7 @@ private:
if(!is_positive(traits_.compute_scalar_product_3_object()(old_n, new_n))) if(!is_positive(traits_.compute_scalar_product_3_object()(old_n, new_n)))
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG #ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Moving vertex would result in the inversion of a face normal!" << std::endl; std::cout << "Moving vertex would result in the inversion of a face normal!" << std::endl;
#endif #endif
return true; return true;
@ -655,7 +707,7 @@ private:
const Point& new_pos) const const Point& new_pos) const
{ {
// check if the minimum angle of the star has not deteriorated // check if the minimum angle of the star has not deteriorated
double old_min_angle = CGAL_PI; FT old_min_angle = CGAL_PI;
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_)) for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{ {
const Point_ref old_pos = get(vpmap_, v); const Point_ref old_pos = get(vpmap_, v);
@ -683,7 +735,7 @@ private:
get_radian_angle(Vector(lpt, rpt), Vector(lpt, new_pos), traits_) < old_min_angle || get_radian_angle(Vector(lpt, rpt), Vector(lpt, new_pos), traits_) < old_min_angle ||
get_radian_angle(Vector(rpt, new_pos), Vector(rpt, lpt), traits_) < old_min_angle) get_radian_angle(Vector(rpt, new_pos), Vector(rpt, lpt), traits_) < old_min_angle)
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG #ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
const Point_ref old_pos = get(vpmap_, v); const Point_ref old_pos = get(vpmap_, v);
std::cout << "deterioration of min angle in the star!" << std::endl; std::cout << "deterioration of min angle in the star!" << std::endl;

6
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_evaluation.h
View File

@ -61,7 +61,7 @@ public:
angles_.push_back(traits_.compute_approximate_angle_3_object()(a, b, c)); angles_.push_back(traits_.compute_approximate_angle_3_object()(a, b, c));
} }
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "angles_ size = " << angles_.size() << std::endl; std::cout << "angles_ size = " << angles_.size() << std::endl;
#endif #endif
} }
@ -79,7 +79,7 @@ public:
for(face_descriptor f : faces(mesh_)) for(face_descriptor f : faces(mesh_))
areas_.push_back(face_area(f, mesh_)); areas_.push_back(face_area(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "areas_ size = " << areas_.size() << std::endl; std::cout << "areas_ size = " << areas_.size() << std::endl;
#endif #endif
} }
@ -97,7 +97,7 @@ public:
for(face_descriptor f : faces(mesh_)) for(face_descriptor f : faces(mesh_))
aspect_ratios_.push_back(CGAL::Polygon_mesh_processing::face_aspect_ratio(f, mesh_)); aspect_ratios_.push_back(CGAL::Polygon_mesh_processing::face_aspect_ratio(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "aspect_ratios_ size = " << aspect_ratios_.size() << std::endl; std::cout << "aspect_ratios_ size = " << aspect_ratios_.size() << std::endl;
#endif #endif
} }

709
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/remove_degeneracies.h
View File

@ -1,709 +0,0 @@
// Copyright (c) 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 <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/boost/graph/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>
#include <iostream>
#include <set>
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
#include <sstream>
#include <fstream>
#endif
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template <typename TriangleMesh, typename VPM, typename ECM, typename Traits>
std::array<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor, 2>
is_badly_shaped(const typename boost::graph_traits<TriangleMesh>::face_descriptor f,
TriangleMesh& tmesh,
const VPM& vpm,
const ECM& ecm,
const Traits& gt,
const double cap_threshold, // angle over 160° ==> cap
const double needle_threshold, // longest edge / shortest edge over this ratio ==> needle
const double collapse_length_threshold) // max length of edges allowed to be collapsed
{
namespace PMP = CGAL::Polygon_mesh_processing;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
const halfedge_descriptor null_h = boost::graph_traits<TriangleMesh>::null_halfedge();
halfedge_descriptor res = PMP::is_needle_triangle_face(f, tmesh, needle_threshold,
parameters::vertex_point_map(vpm)
.geom_traits(gt));
if(res != null_h && !get(ecm, edge(res, tmesh)))
{
// don't want to collapse edges that are too large
if(collapse_length_threshold == 0 ||
edge_length(res, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)) <= collapse_length_threshold)
{
return make_array(res, null_h);
}
}
else // let's not make it possible to have a face be both a cap and a needle (for now)
{
res = PMP::is_cap_triangle_face(f, tmesh, cap_threshold, parameters::vertex_point_map(vpm).geom_traits(gt));
if(res != null_h && !get(ecm, edge(res, tmesh)))
return make_array(null_h, res);
}
return make_array(null_h, null_h);
}
template <typename TriangleMesh, typename EdgeContainer,
typename VPM, typename ECM, typename Traits>
void collect_badly_shaped_triangles(const typename boost::graph_traits<TriangleMesh>::face_descriptor f,
TriangleMesh& tmesh,
const VPM& vpm,
const ECM& ecm,
const Traits& gt,
const double cap_threshold, // angle over this threshold (as a cosine) ==> cap
const double needle_threshold, // longest edge / shortest edge over this ratio ==> needle
const double collapse_length_threshold, // max length of edges allowed to be collapsed
EdgeContainer& edges_to_collapse,
EdgeContainer& edges_to_flip)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
std::array<halfedge_descriptor, 2> res = is_badly_shaped(f, tmesh, vpm, ecm, gt, cap_threshold,
needle_threshold, collapse_length_threshold);
if(res[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "add new needle: " << edge(res[0], tmesh) << std::endl;
#endif
edges_to_collapse.insert(edge(res[0], tmesh));
}
else // let's not make it possible to have a face be both a cap and a needle (for now)
{
if(res[1] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "add new cap: " << edge(res[1],tmesh) << std::endl;
#endif
edges_to_flip.insert(edge(res[1], tmesh));
}
}
}
/*
// Following Ronfard et al. 96 we look at variation of the normal after the collapse
// the collapse must be topologically valid
template <class TriangleMesh, class NamedParameters>
bool is_collapse_geometrically_valid(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h,
const TriangleMesh& tmesh,
const NamedParameters& np)
{
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VPM;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tmesh));
Traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
/// @todo handle boundary edges
h = opposite(h, tmesh); // Euler::collapse edge keeps the target and removes the source
// source is kept, target is removed
CGAL_assertion(target(h, tmesh) == vertex_removed);
Point_ref kept = get(vpm, source(h, tmesh));
Point_ref removed= get(vpm, target(h, tmesh));
// consider triangles incident to the vertex removed
halfedge_descriptor stop = prev(opposite(h, tmesh), tmesh);
halfedge_descriptor hi = opposite(next(h, tmesh), tmesh);
std::vector<halfedge_descriptor> triangles;
while(hi != stop)
{
if(!is_border(hi, tmesh))
{
Point_ref a = get(vpm, target(next(hi, tmesh), tmesh));
Point_ref b = get(vpm, source(hi, tmesh));
//ack a-b-point_remove and a-b-point_kept has a compatible orientation
/// @todo use a predicate
typename Traits::Vector_3 n1 = gt.construct_cross_product_vector_3_object()(removed-a, b-a);
typename Traits::Vector_3 n2 = gt.construct_cross_product_vector_3_object()(kept-a, b-a);
if(gt.compute_scalar_product_3_object()(n1, n2) <= 0)
return false;
}
hi = opposite(next(hi, tmesh), tmesh);
}
return true;
}
*/
template <class TriangleMesh, typename VPM, typename Traits>
boost::optional<double> get_collapse_volume(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h,
const TriangleMesh& tmesh,
const VPM& vpm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename Traits::Vector_3 Vector_3;
const typename Traits::Point_3 origin(ORIGIN);
/// @todo handle boundary edges
h = opposite(h, tmesh); // Euler::collapse edge keeps the target and removes the source
// source is kept, target is removed
Point_ref kept = get(vpm, source(h, tmesh));
Point_ref removed= get(vpm, target(h, tmesh));
// init volume with incident triangles (reversed orientation
double delta_vol = volume(removed, kept, get(vpm, target(next(h, tmesh), tmesh)), origin) +
volume(kept, removed, get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh)), origin);
// consider triangles incident to the vertex removed
halfedge_descriptor stop = prev(opposite(h, tmesh), tmesh);
halfedge_descriptor hi = opposite(next(h, tmesh), tmesh);
std::vector<halfedge_descriptor> triangles;
while(hi != stop)
{
if(!is_border(hi, tmesh))
{
Point_ref a = get(vpm, target(next(hi, tmesh), tmesh));
Point_ref b = get(vpm, source(hi, tmesh));
//ack a-b-point_remove and a-b-point_kept has a compatible orientation
/// @todo use a predicate
Vector_3 v_ab = gt.construct_vector_3_object()(a, b);
Vector_3 v_ar = gt.construct_vector_3_object()(a, removed);
Vector_3 v_ak = gt.construct_vector_3_object()(a, kept);
Vector_3 n1 = gt.construct_cross_product_vector_3_object()(v_ar, v_ab);
Vector_3 n2 = gt.construct_cross_product_vector_3_object()(v_ak, v_ab);
if(gt.compute_scalar_product_3_object()(n1, n2) <= 0)
return boost::none;
delta_vol += volume(b, a, removed, origin) + volume(a, b, kept, origin); // opposite orientation
}
hi = opposite(next(hi, tmesh), tmesh);
}
return CGAL::abs(delta_vol);
}
template <typename TriangleMesh, typename VPM, typename VCM, typename Traits>
typename boost::graph_traits<TriangleMesh>::halfedge_descriptor
get_best_edge_orientation(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
const TriangleMesh& tmesh,
const VPM& vpm,
const VCM& vcm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
halfedge_descriptor h = halfedge(e, tmesh), ho = opposite(h, tmesh);
CGAL_assertion(!get(vcm, source(h, tmesh)) || !get(vcm, target(h, tmesh)));
boost::optional<double> dv1 = get_collapse_volume(h, tmesh, vpm, gt);
boost::optional<double> dv2 = get_collapse_volume(ho, tmesh, vpm, gt);
// the resulting point of the collapse of a halfedge is the target of the halfedge before collapse
if(get(vcm, source(h, tmesh)))
return dv2 != boost::none ? ho
: boost::graph_traits<TriangleMesh>::null_halfedge();
if(get(vcm, target(h, tmesh)))
return dv1 != boost::none ? h
: boost::graph_traits<TriangleMesh>::null_halfedge();
if(dv1 != boost::none)
{
if(dv2 != boost::none)
return (*dv1 < *dv2) ? h : ho;
return h;
}
if(dv2 != boost::none)
return ho;
return boost::graph_traits<TriangleMesh>::null_halfedge();
}
// adapted from triangulate_faces
template <typename TriangleMesh, typename VPM, typename Traits>
bool should_flip(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
const TriangleMesh& tmesh,
const VPM& vpm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename Traits::Vector_3 Vector_3;
CGAL_precondition(!is_border(e, tmesh));
halfedge_descriptor h = halfedge(e, tmesh);
Point_ref p0 = get(vpm, target(h, tmesh));
Point_ref p1 = get(vpm, target(next(h, tmesh), tmesh));
Point_ref p2 = get(vpm, source(h, tmesh));
Point_ref p3 = get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh));
/* Chooses the diagonal that will split the quad in two triangles that maximize
* the scalar product of of the un-normalized normals of the two triangles.
* The lengths of the un-normalized normals (computed using cross-products of two vectors)
* are proportional to the area of the triangles.
* Maximize the scalar product of the two normals will avoid skinny triangles,
* and will also taken into account the cosine of the angle between the two normals.
* In particular, if the two triangles are oriented in different directions,
* the scalar product will be negative.
*/
// CGAL::cross_product(p2-p1, p3-p2) * CGAL::cross_product(p0-p3, p1-p0);
// CGAL::cross_product(p1-p0, p1-p2) * CGAL::cross_product(p3-p2, p3-p0);
const Vector_3 v01 = gt.construct_vector_3_object()(p0, p1);
const Vector_3 v12 = gt.construct_vector_3_object()(p1, p2);
const Vector_3 v23 = gt.construct_vector_3_object()(p2, p3);
const Vector_3 v30 = gt.construct_vector_3_object()(p3, p0);
const double p1p3 = gt.compute_scalar_product_3_object()(
gt.construct_cross_product_vector_3_object()(v12, v23),
gt.construct_cross_product_vector_3_object()(v30, v01));
const Vector_3 v21 = gt.construct_opposite_vector_3_object()(v12);
const Vector_3 v03 = gt.construct_opposite_vector_3_object()(v30);
const double p0p2 = gt.compute_scalar_product_3_object()(
gt.construct_cross_product_vector_3_object()(v01, v21),
gt.construct_cross_product_vector_3_object()(v23, v03));
return p0p2 <= p1p3;
}
} // namespace internal
namespace experimental {
// @todo check what to use as priority queue with removable elements, set might not be optimal
template <typename FaceRange, typename TriangleMesh, typename NamedParameters>
bool remove_almost_degenerate_faces(const FaceRange& face_range,
TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold,
const NamedParameters& np)
{
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef Constant_property_map<vertex_descriptor, bool> Default_VCM;
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
NamedParameters,
Default_VCM>::type VCM;
VCM vcm_np = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), Default_VCM(false));
typedef Constant_property_map<edge_descriptor, bool> Default_ECM;
typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
NamedParameters,
Default_ECM>::type ECM;
ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), Default_ECM(false));
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VPM;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tmesh));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
Traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
// Vertex property map that combines the VCM and the fact that extremities of a constrained edge should be constrained
typedef CGAL::dynamic_vertex_property_t<bool> Vertex_property_tag;
typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type DVCM;
DVCM vcm = get(Vertex_property_tag(), tmesh);
for(face_descriptor f : face_range)
{
if(f == boost::graph_traits<TriangleMesh>::null_face())
continue;
for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh))
{
if(get(ecm, edge(h, tmesh)))
{
put(vcm, source(h, tmesh), true);
put(vcm, target(h, tmesh), true);
}
else if(get(vcm_np, target(h, tmesh)))
{
put(vcm, target(h, tmesh), true);
}
}
}
// Start the process of removing bad elements
std::set<edge_descriptor> edges_to_collapse;
std::set<edge_descriptor> edges_to_flip;
// @todo could probably do something a bit better by looping edges, consider the incident faces
// f1 / f2 and look at f1 if f1<f2, and the edge is smaller than the two other edges...
for(face_descriptor f : face_range)
internal::collect_badly_shaped_triangles(f, tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
int iter = 0;
#endif
for(;;)
{
bool something_was_done = false;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << edges_to_collapse.size() << " needles and " << edges_to_flip.size() << " caps" << std::endl;
std::ostringstream oss;
oss << "degen_cleaning_iter_" << iter++ << ".off";
std::ofstream out(oss.str().c_str());
out << std::setprecision(17);
out << tmesh;
out.close();
#endif
if(edges_to_collapse.empty() && edges_to_flip.empty())
return true;
/// @todo maybe using a priority queue handling the more almost degenerate elements should be used
std::set<edge_descriptor> next_edges_to_collapse;
std::set<edge_descriptor> next_edges_to_flip;
// treat needles
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
int kk=0;
std::ofstream(std::string("tmp/n-00000.off")) << tmesh;
#endif
while(!edges_to_collapse.empty())
{
edge_descriptor e = *edges_to_collapse.begin();
edges_to_collapse.erase(edges_to_collapse.begin());
CGAL_assertion(!get(ecm, e));
if(get(vcm, source(e, tmesh)) && get(vcm, target(e, tmesh)))
continue;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << " treat needle: " << e << " (" << tmesh.point(source (e, tmesh))
<< " --- " << tmesh.point(target(e, tmesh)) << ")" << std::endl;
#endif
if(CGAL::Euler::does_satisfy_link_condition(e, tmesh))
{
// the following edges are removed by the collapse
halfedge_descriptor h = halfedge(e, tmesh);
CGAL_assertion(!is_border(h, tmesh)); // because extracted from a face
std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[0] != h)
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: Needle criteria no longer verified " << tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
// the opposite edge might also have been inserted in the set and might still be a needle
h = opposite(h, tmesh);
if(is_border(h, tmesh))
continue;
nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold,
collapse_length_threshold);
if(nc[0] != h)
continue;
}
for(int i=0; i<2; ++i)
{
if(!is_border(h, tmesh))
{
edge_descriptor pe = edge(prev(h, tmesh), tmesh);
edges_to_flip.erase(pe);
next_edges_to_collapse.erase(pe);
edges_to_collapse.erase(pe);
}
h = opposite(h, tmesh);
}
// pick the orientation of edge to keep the vertex minimizing the volume variation
h = internal::get_best_edge_orientation(e, tmesh, vpm, vcm, gt);
if(h == boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: geometrically invalid edge collapse! "
<< tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
next_edges_to_collapse.insert(e);
continue;
}
edges_to_flip.erase(e);
next_edges_to_collapse.erase(e); // for edges added in faces incident to a vertex kept after a collapse
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cerr << " " << kk << " -- Collapsing " << tmesh.point(source(h, tmesh)) << " "
<< tmesh.point(target(h, tmesh)) << std::endl;
#endif
// moving to the midpoint is not a good idea. On a circle for example you might endpoint with
// a bad geometry because you iteratively move one point
// auto mp = midpoint(tmesh.point(source(h, tmesh)), tmesh.point(target(h, tmesh)));
vertex_descriptor v = Euler::collapse_edge(edge(h, tmesh), tmesh);
//tmesh.point(v) = mp;
// examine all faces incident to the vertex kept
for(halfedge_descriptor hv : halfedges_around_target(v, tmesh))
{
if(!is_border(hv, tmesh))
{
internal::collect_badly_shaped_triangles(face(hv, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
}
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::string nb = std::to_string(++kk);
if(kk<10) nb = std::string("0")+nb;
if(kk<100) nb = std::string("0")+nb;
if(kk<1000) nb = std::string("0")+nb;
if(kk<10000) nb = std::string("0")+nb;
std::ofstream(std::string("tmp/n-")+nb+std::string(".off")) << tmesh;
#endif
something_was_done = true;
}
else
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: uncollapsable edge! " << tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
next_edges_to_collapse.insert(e);
}
}
// treat caps
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
kk=0;
std::ofstream(std::string("tmp/c-000.off")) << tmesh;
#endif
while(!edges_to_flip.empty())
{
edge_descriptor e = *edges_to_flip.begin();
edges_to_flip.erase(edges_to_flip.begin());
CGAL_assertion(!get(ecm, e));
if(get(vcm, source(e, tmesh)) && get(vcm, target(e, tmesh)))
continue;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "treat cap: " << e << " (" << tmesh.point(source(e, tmesh))
<< " --- " << tmesh.point(target(e, tmesh)) << ")" << std::endl;
#endif
halfedge_descriptor h = halfedge(e, tmesh);
std::array<halfedge_descriptor,2> nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold,
collapse_length_threshold);
// First check the triangle is still a cap
if(nc[1] != h)
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: Cap criteria no longer verified " << tmesh.point(source(e, tmesh)) << " --- "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
// the opposite edge might also have been inserted in the set and might still be a cap
h = opposite(h, tmesh);
if(is_border(h, tmesh))
continue;
nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[1] != h)
continue;
}
// special case on the border
if(is_border(opposite(h, tmesh), tmesh))
{
// remove the triangle
edges_to_flip.erase(edge(prev(h, tmesh), tmesh));
edges_to_flip.erase(edge(next(h, tmesh), tmesh));
next_edges_to_collapse.erase(edge(prev(h, tmesh), tmesh));
next_edges_to_collapse.erase(edge(next(h, tmesh), tmesh));
Euler::remove_face(h, tmesh);
something_was_done = true;
continue;
}
// condition for the flip to be valid (the edge to be created does not already exist)
if(!halfedge(target(next(h, tmesh), tmesh),
target(next(opposite(h, tmesh), tmesh), tmesh), tmesh).second)
{
if(!internal::should_flip(e, tmesh, vpm, gt))
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "Flipping prevented: not the best diagonal" << std::endl;
#endif
next_edges_to_flip.insert(e);
continue;
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "Flipping" << std::endl;
#endif
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cerr << "step " << kk << "\n";
std::cerr << " Flipping " << tmesh.point(source(h, tmesh)) << " "
<< tmesh.point(target(h, tmesh)) << std::endl;
#endif
Euler::flip_edge(h, tmesh);
CGAL_assertion(edge(h, tmesh) == e);
// handle face updates
for(int i=0; i<2; ++i)
{
CGAL_assertion(!is_border(h, tmesh));
std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[1] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
if(edge(nc[1], tmesh) != e)
next_edges_to_flip.insert(edge(nc[1], tmesh));
}
else
{
if(nc[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
next_edges_to_collapse.insert(edge(nc[0], tmesh));
}
}
h = opposite(h, tmesh);
}
something_was_done = true;
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
else
{
std::cerr << "Warning: unflippable edge! " << tmesh.point(source(h, tmesh)) << " --- "
<< tmesh.point(target(h, tmesh)) << std::endl;
next_edges_to_flip.insert(e);
}
#endif
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::string nb = std::to_string(++kk);
if(kk<10) nb = std::string("0")+nb;
if(kk<100) nb = std::string("0")+nb;
if(kk<1000) nb = std::string("0")+nb;
if(kk<10000) nb = std::string("0")+nb;
std::ofstream(std::string("tmp/c-")+nb+std::string(".off")) << tmesh;
#endif
}
std::swap(edges_to_collapse, next_edges_to_collapse);
std::swap(edges_to_flip, next_edges_to_flip);
if(!something_was_done)
return false;
}
return false;
}
template <typename FaceRange, typename TriangleMesh>
bool remove_almost_degenerate_faces(const FaceRange& face_range,
TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold)
{
return remove_almost_degenerate_faces(face_range, tmesh,
cap_threshold, needle_threshold, collapse_length_threshold,
CGAL::parameters::all_default());
}
template <typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
bool remove_almost_degenerate_faces(TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold,
const CGAL_PMP_NP_CLASS& np)
{
return remove_almost_degenerate_faces(faces(tmesh), tmesh, cap_threshold, needle_threshold,
collapse_length_threshold, np);
}
template<class TriangleMesh>
bool remove_almost_degenerate_faces(TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold)
{
return remove_almost_degenerate_faces(faces(tmesh), tmesh,
cap_threshold, needle_threshold, collapse_length_threshold,
CGAL::parameters::all_default());
}
} // namespace experimental
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H

461
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/manifoldness.h Normal file
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

94
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h
View File

@ -17,23 +17,99 @@
#include <CGAL/disable_warnings.h> #include <CGAL/disable_warnings.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/property_map.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h> #include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/algorithm.h>
#include <set>
#include <boost/dynamic_bitset.hpp>
#include <CGAL/algorithm.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/property_map.h>
#include <boost/dynamic_bitset.hpp>
#include <boost/range/size.hpp> #include <boost/range/size.hpp>
#include <boost/range/value_type.hpp> #include <boost/range/value_type.hpp>
#include <boost/range/reference.hpp> #include <boost/range/reference.hpp>
namespace CGAL #include <set>
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
// \ingroup PMP_repairing_grp
//
// Adds the vertices and faces of a mesh into a (possibly non-empty) soup.
//
// @tparam PolygonMesh a model of `FaceListGraph` with an internal point property map
// @tparam PointRange a model of the concepts `RandomAccessContainer` and
// `BackInsertionSequence` whose value type is the point type
// @tparam PolygonRange a model of the concepts `RandomAccessContainer` and `BackInsertionSequence` whose
// `value_type` is itself a model of the concepts `RandomAccessContainer` and
// `BackInsertionSequence` whose `value_type` is `std::size_t`.
//
// @param mesh the mesh whose faces are being put in the soup
// @param points points of the soup of polygons
// @param polygons each element in the vector describes a polygon using the index of the points in `points`
//
// \sa `CGAL::Polygon_mesh_processing::orient_polygon_soup()`
// \sa `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`
// \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
//
template<typename PolygonMesh, typename PointRange, typename PolygonRange, typename NamedParameters>
void polygon_mesh_to_polygon_soup(const PolygonMesh& mesh,
PointRange& points,
PolygonRange& polygons,
const NamedParameters& np)
{ {
namespace Polygon_mesh_processing typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
{ typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
namespace internal typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
using parameters::choose_parameter;
using parameters::get_parameter;
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPM;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, mesh));
typedef CGAL::dynamic_vertex_property_t<std::size_t> Vertex_index;
typedef typename boost::property_map<PolygonMesh, Vertex_index>::const_type VIM;
VIM vim = get(Vertex_index(), mesh);
typedef typename boost::range_value<PolygonRange>::type Polygon;
std::size_t index = points.size(); // so that multiple meshes can be put into the same soup
points.reserve(points.size() + vertices(mesh).size());
polygons.reserve(polygons.size() + faces(mesh).size());
for(const vertex_descriptor v : vertices(mesh))
{
points.push_back(get(vpm, v));
put(vim, v, index++);
}
for(const face_descriptor f : faces(mesh))
{
Polygon polygon;
for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, mesh), mesh))
polygon.push_back(get(vim, target(h, mesh)));
polygons.push_back(polygon);
}
}
template<typename PolygonMesh, typename PointRange, typename PolygonRange>
void polygon_mesh_to_polygon_soup(const PolygonMesh& mesh,
PointRange& points,
PolygonRange& polygons)
{ {
return polygon_mesh_to_polygon_soup(mesh, points, polygons, CGAL::parameters::all_default());
}
// -------------------------------------------------------------------------------------------------
template <typename PM template <typename PM
, typename PointRange , typename PointRange
, typename PolygonRange> , typename PolygonRange>

2397
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_degeneracies.h Normal file
View File

File diff suppressed because it is too large Load Diff

1865
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remove_self_intersections.h Normal file
View File

File diff suppressed because it is too large Load Diff

2938
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h
View File

File diff suppressed because it is too large Load Diff

153
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_predicates.h
View File

@ -23,6 +23,7 @@
#include <CGAL/boost/graph/iterator.h> #include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/helpers.h> #include <CGAL/boost/graph/helpers.h>
#include <boost/range/has_range_iterator.hpp>
#include <boost/graph/graph_traits.hpp> #include <boost/graph/graph_traits.hpp>
#include <limits> #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()); 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 /// \ingroup PMP_repairing_grp
/// checks whether a triangle face is degenerate. /// checks whether a triangle face is degenerate.
/// A triangle face is considered degenerate if the geometric positions of its vertices are collinear. /// 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()); 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 /// \ingroup PMP_repairing_grp
/// checks whether a triangle face is needle. /// 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. /// A triangle is said to be a <i>needle</i> if its longest edge is much longer than its shortest edge.

20
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h
View File

@ -173,7 +173,7 @@ void smooth_mesh(const FaceRange& faces,
const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true); const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true);
VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
get(Vertex_property_tag(), tmesh)); get(Vertex_property_tag(), tmesh));
// If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization // If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization
if((std::is_same<VCMap, Default_VCMap>::value)) if((std::is_same<VCMap, Default_VCMap>::value))
@ -183,7 +183,7 @@ void smooth_mesh(const FaceRange& faces,
} }
ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
Constant_property_map<edge_descriptor, bool>(false)); Constant_property_map<edge_descriptor, bool>(false));
// a constrained edge has constrained extremities // a constrained edge has constrained extremities
for(face_descriptor f : faces) for(face_descriptor f : faces)
@ -235,8 +235,16 @@ void smooth_mesh(const FaceRange& faces,
for(unsigned int i=0; i<nb_iterations; ++i) for(unsigned int i=0; i<nb_iterations; ++i)
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Iteration #" << i << std::endl;
#endif
if(use_area_smoothing) if(use_area_smoothing)
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Smooth areas..." << std::endl;
#endif
// First apply area smoothing... // First apply area smoothing...
area_smoother.optimize(use_safety_constraints /*check for bad faces*/, area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
false /*apply moves as soon as they're calculated*/, false /*apply moves as soon as they're calculated*/,
@ -245,7 +253,7 @@ void smooth_mesh(const FaceRange& faces,
{ {
if(use_safety_constraints && does_self_intersect(tmesh)) if(use_safety_constraints && does_self_intersect(tmesh))
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n"; std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
#endif #endif
break; break;
@ -261,6 +269,10 @@ void smooth_mesh(const FaceRange& faces,
// ... then angle smoothing // ... then angle smoothing
if(use_angle_smoothing) if(use_angle_smoothing)
{ {
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Smooth angles..." << std::endl;
#endif
angle_smoother.optimize(use_safety_constraints /*check for bad faces*/, angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
true /*apply all moves at once*/, true /*apply all moves at once*/,
use_safety_constraints /*check if the min angle is improved*/); use_safety_constraints /*check if the min angle is improved*/);
@ -269,7 +281,7 @@ void smooth_mesh(const FaceRange& faces,
{ {
if(use_safety_constraints && does_self_intersect(tmesh)) if(use_safety_constraints && does_self_intersect(tmesh))
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n"; std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
#endif #endif
break; break;

4
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h
View File

@ -149,7 +149,7 @@ void smooth_shape(const FaceRange& faces,
for(unsigned int iter=0; iter<nb_iterations; ++iter) for(unsigned int iter=0; iter<nb_iterations; ++iter)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "iteration #" << iter << std::endl; std::cout << "iteration #" << iter << std::endl;
#endif #endif
@ -161,7 +161,7 @@ void smooth_shape(const FaceRange& faces,
} }
else else
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Failed to solve system!" << std::endl; std::cerr << "Failed to solve system!" << std::endl;
#endif #endif
break; break;

10
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h
View File

@ -843,7 +843,7 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
template <typename PolygonMesh, template <typename PolygonMesh,
typename HalfedgePairsRange> typename HalfedgePairsRange>
std::size_t stitch_borders(PolygonMesh& pmesh, std::size_t stitch_borders(PolygonMesh& pmesh,
const HalfedgePairsRange& hedge_pairs_to_stitch) const HalfedgePairsRange& hedge_pairs_to_stitch)
{ {
return stitch_borders(pmesh, hedge_pairs_to_stitch, CGAL::parameters::all_default()); return stitch_borders(pmesh, hedge_pairs_to_stitch, CGAL::parameters::all_default());
} }
@ -890,9 +890,9 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
halfedge_descriptor; halfedge_descriptor;
std::vector< std::pair<halfedge_descriptor, halfedge_descriptor> > hedge_pairs_to_stitch; std::vector< std::pair<halfedge_descriptor, halfedge_descriptor> > hedge_pairs_to_stitch;
typedef typename GetVertexPointMap<PolygonMesh, CGAL_PMP_NP_CLASS>::const_type VPMap; typedef typename GetVertexPointMap<PolygonMesh, CGAL_PMP_NP_CLASS>::const_type VPM;
VPMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh)); get_const_property_map(vertex_point, pmesh));
#ifdef CGAL_PMP_STITCHING_DEBUG #ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "------- Stitch cycles..." << std::endl; std::cout << "------- Stitch cycles..." << std::endl;
@ -907,7 +907,7 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
internal::collect_duplicated_stitchable_boundary_edges(pmesh, internal::collect_duplicated_stitchable_boundary_edges(pmesh,
std::back_inserter(hedge_pairs_to_stitch), std::back_inserter(hedge_pairs_to_stitch),
internal::Less_for_halfedge<PolygonMesh, VPMap>(pmesh, vpm), internal::Less_for_halfedge<PolygonMesh, VPM>(pmesh, vpm),
vpm, np); vpm, np);
res += stitch_borders(pmesh, hedge_pairs_to_stitch, np); res += stitch_borders(pmesh, hedge_pairs_to_stitch, np);

4
Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h
View File

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

3
Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
View File

@ -98,6 +98,7 @@ endif()
create_single_source_cgal_program("test_pmp_manifoldness.cpp") create_single_source_cgal_program("test_pmp_manifoldness.cpp")
create_single_source_cgal_program("test_mesh_smoothing.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_remove_caps_needles.cpp")
create_single_source_cgal_program("test_pmp_remove_self_intersections.cpp")
if( TBB_FOUND ) if( TBB_FOUND )
CGAL_target_use_TBB(test_pmp_distance) CGAL_target_use_TBB(test_pmp_distance)
@ -115,6 +116,8 @@ find_package(Ceres QUIET)
if(TARGET ceres) if(TARGET ceres)
target_compile_definitions( test_mesh_smoothing PRIVATE CGAL_PMP_USE_CERES_SOLVER ) target_compile_definitions( test_mesh_smoothing PRIVATE CGAL_PMP_USE_CERES_SOLVER )
target_link_libraries( test_mesh_smoothing PRIVATE ceres ) 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) endif(TARGET ceres)
if(BUILD_TESTING) if(BUILD_TESTING)

131
Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_remove_self_intersections.cpp Normal file
View File

@ -0,0 +1,131 @@
// #define CGAL_PMP_SMOOTHING_DEBUG
#define CGAL_PMP_COMPUTE_NORMAL_DEBUG
#define CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
#define CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
#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;
namespace CP = ::CGAL::parameters;
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* mesh_filename,
const char* mesh_selection_filename = nullptr)
{
std::cout << std::endl << "---------------" << std::endl;
std::cout << "Test " << mesh_filename << std::endl;
if(mesh_selection_filename)
std::cout << "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::list<face_descriptor> selected_faces;
if(mesh_selection_filename)
{
std::ifstream selection_input(mesh_selection_filename);
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::experimental::remove_self_intersections(selected_faces, mesh);
}
else
{
PMP::experimental::remove_self_intersections(mesh);
}
std::ofstream out("post_repair.off");
out.precision(17);
out << mesh;
out.close();
assert(CGAL::is_valid_polygon_mesh(mesh));
assert(!PMP::does_self_intersect(mesh));
}
void fix_local_self_intersections(const char* fname)
{
std::cout << std::endl << "-----LOCAL------" << 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);
}
PMP::experimental::remove_self_intersections(mesh, CP::apply_per_connected_component(true));
std::ofstream out("post_local_repair.off");
out.precision(17);
out << mesh;
out.close();
assert(CGAL::is_valid_polygon_mesh(mesh));
assert(!PMP::does_self_intersect(mesh));
}
int main()
{
std::cout.precision(17);
std::cerr.precision(17);
fix_local_self_intersections("data_repair/brain.off");
#if 0
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");
// fix_self_intersections("data_repair/hemispheres.off", "data_repair/hemispheres-half.selection.txt");
return 0;
}

2
Polygon_mesh_processing/test/Polygon_mesh_processing/test_remove_caps_needles.cpp
View File

@ -3,7 +3,7 @@
#include <CGAL/Polygon_mesh_processing/self_intersections.h> #include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <fstream> #include <fstream>
#include <CGAL/Polygon_mesh_processing/internal/remove_degeneracies.h> #include <CGAL/Polygon_mesh_processing/remove_degeneracies.h>
#include <iostream> #include <iostream>
#include <vector> #include <vector>

18
Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp
View File

@ -1,4 +1,4 @@
#define CGAL_PMP_SMOOTHING_VERBOSE #define CGAL_PMP_SMOOTHING_DEBUG
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@ -31,7 +31,7 @@ bool equal_doubles(double d1, double d2, double e)
template <typename Mesh> template <typename Mesh>
void test_implicit_constrained_devil(Mesh mesh) void test_implicit_constrained_devil(Mesh mesh)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "-- test_implicit_constrained_devil --" << std::endl; std::cout << "-- test_implicit_constrained_devil --" << std::endl;
#endif #endif
@ -69,7 +69,7 @@ void test_implicit_constrained_devil(Mesh mesh)
++i; ++i;
} }
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::ofstream out("output_implicit_constrained_devil.off"); std::ofstream out("output_implicit_constrained_devil.off");
out << mesh; out << mesh;
out.close(); out.close();
@ -79,7 +79,7 @@ void test_implicit_constrained_devil(Mesh mesh)
template <typename Mesh> template <typename Mesh>
void test_implicit_constrained_elephant(Mesh mesh) void test_implicit_constrained_elephant(Mesh mesh)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "-- test_implicit_constrained_elephant --" << std::endl; std::cout << "-- test_implicit_constrained_elephant --" << std::endl;
#endif #endif
@ -117,7 +117,7 @@ void test_implicit_constrained_elephant(Mesh mesh)
++i; ++i;
} }
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::ofstream out("output_implicit_constrained_elephant.off"); std::ofstream out("output_implicit_constrained_elephant.off");
out << mesh; out << mesh;
out.close(); out.close();
@ -127,14 +127,14 @@ void test_implicit_constrained_elephant(Mesh mesh)
template <typename Mesh> template <typename Mesh>
void test_curvature_flow_time_step(Mesh mesh) void test_curvature_flow_time_step(Mesh mesh)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "-- test_curvature_flow_time_step --" << std::endl; std::cout << "-- test_curvature_flow_time_step --" << std::endl;
#endif #endif
const double time_step = 1e-15; const double time_step = 1e-15;
PMP::smooth_shape(mesh, time_step); PMP::smooth_shape(mesh, time_step);
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::ofstream out("output_devil_time_step.off"); std::ofstream out("output_devil_time_step.off");
out << mesh; out << mesh;
out.close(); out.close();
@ -144,14 +144,14 @@ void test_curvature_flow_time_step(Mesh mesh)
template <typename Mesh> template <typename Mesh>
void test_curvature_flow(Mesh mesh) void test_curvature_flow(Mesh mesh)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "-- test_curvature_flow --" << std::endl; std::cout << "-- test_curvature_flow --" << std::endl;
#endif #endif
const double time_step = 1.0; const double time_step = 1.0;
PMP::smooth_shape(mesh, time_step); PMP::smooth_shape(mesh, time_step);
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_DEBUG
std::ofstream out("output_precision_elephant.off"); std::ofstream out("output_precision_elephant.off");
out << mesh; out << mesh;
out.close(); out.close();

4
Polyhedron/demo/Polyhedron/Plugins/PMP/Repair_polyhedron_plugin.cpp
View File

@ -16,7 +16,7 @@
#include <CGAL/Polygon_mesh_processing/repair.h> #include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/internal/repair_extra.h> #include <CGAL/Polygon_mesh_processing/internal/repair_extra.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h> #include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/internal/remove_degeneracies.h> #include <CGAL/Polygon_mesh_processing/remove_degeneracies.h>
#include "ui_RemoveNeedlesDialog.h" #include "ui_RemoveNeedlesDialog.h"
@ -209,7 +209,7 @@ void Polyhedron_demo_repair_polyhedron_plugin::on_actionRemoveSelfIntersections_
if (poly_item) if (poly_item)
{ {
bool solved = bool solved =
CGAL::Polygon_mesh_processing::remove_self_intersections( CGAL::Polygon_mesh_processing::experimental::remove_self_intersections(
*poly_item->polyhedron()); *poly_item->polyhedron());
if (!solved) if (!solved)
CGAL::Three::Three::information(tr("Some self-intersection could not be fixed")); CGAL::Three::Three::information(tr("Some self-intersection could not be fixed"));

28
Polyhedron/demo/Polyhedron/Scene_polygon_soup_item.cpp
View File

@ -43,6 +43,7 @@
#include <boost/accumulators/statistics/max.hpp> #include <boost/accumulators/statistics/max.hpp>
#include <boost/accumulators/statistics/median.hpp> #include <boost/accumulators/statistics/median.hpp>
#include <iostream>
#include <map> #include <map>
#include <streambuf> #include <streambuf>
@ -366,37 +367,12 @@ void Scene_polygon_soup_item::init_polygon_soup(std::size_t nb_pts, std::size_t
d->oriented = false; d->oriented = false;
} }
#include <iostream>
template<class PolygonMesh> template<class PolygonMesh>
void polygon_mesh_to_soup(PolygonMesh& mesh, Polygon_soup& soup) void polygon_mesh_to_soup(PolygonMesh& mesh, Polygon_soup& soup)
{ {
soup.clear(); soup.clear();
typedef typename boost::property_map<PolygonMesh, boost::vertex_point_t>::type VPMap; CGAL::Polygon_mesh_processing::internal::polygon_mesh_to_polygon_soup(mesh, soup.points, soup.polygons);
VPMap vpmap = get(boost::vertex_point, mesh);
std::map<typename boost::graph_traits<PolygonMesh>::vertex_descriptor, int> vim;
int index=0;
//fill points
for(typename boost::graph_traits<PolygonMesh>::vertex_iterator vit =
vertices(mesh).begin(); vit != vertices(mesh).end(); ++vit)
{
soup.points.push_back(get(vpmap, *vit));
vim.insert(std::make_pair(*vit, index++));
}
//fill triangles
for(typename boost::graph_traits<PolygonMesh>::face_iterator fit =
faces(mesh).begin(); fit != faces(mesh).end(); ++fit)
{
Polygon_soup::Polygon_3 polygon;
for(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor hd :
CGAL::halfedges_around_face(halfedge(*fit, mesh), mesh))
{
polygon.push_back(vim[target(hd, mesh)]);
}
soup.polygons.push_back(polygon);
}
soup.fill_edges(); soup.fill_edges();
} }

1
Polyhedron_IO/include/CGAL/IO/reader_helpers.h
View File

@ -14,6 +14,7 @@
#include <CGAL/array.h> #include <CGAL/array.h>
#include <CGAL/assertions.h> #include <CGAL/assertions.h>
#include <CGAL/Has_member.h> #include <CGAL/Has_member.h>
#include <CGAL/Point_3.h>
#include <boost/mpl/logical.hpp> #include <boost/mpl/logical.hpp>
#include <boost/utility/enable_if.hpp> #include <boost/utility/enable_if.hpp>