Merge pull request #3258 from MaelRL/PMP-Filter_degenerate_polygon_soup_input-GF

PMP: repair polygon soup

BGL/include/CGAL/boost/graph/Seam_mesh.h

@@ -1017,7 +1017,7 @@ public:
   /// of a vertex of the underlying mesh is given by its position
   /// in the container `tm_vds`.
   ///
-  /// \tparam VdContainer must be a model of <a href="http://en.cppreference.com/w/cpp/concept/SequenceContainer"><tt>SequenceContainer</tt></a> (that is, provide
+  /// \tparam VdContainer must be a model of `SequenceContainer` (that is, provide
   ///         the functions: `operator[]` and `at()`).
   ///
   /// \returns one of the halfedges of the seam mesh that is on a seam.
@@ -1062,7 +1062,7 @@ public:
   ///
   /// \returns one of the halfedges of the seam mesh that is on a seam.
   ///
-  /// \tparam VdContainer must be a model of <a href="http://en.cppreference.com/w/cpp/concept/SequenceContainer"><tt>SequenceContainer</tt></a> (that is, provide
+  /// \tparam VdContainer must be a model of `SequenceContainer` (that is, provide
   ///         the functions: `operator[]` and `at()`).
   ///
   /// \pre filename should be the name of a CGAL selection file: edges are

BGL/include/CGAL/boost/graph/parameters_interface.h

@@ -79,6 +79,8 @@ CGAL_add_named_parameter(throw_on_self_intersection_t, throw_on_self_intersectio
 CGAL_add_named_parameter(clip_volume_t, clip_volume, clip_volume)
 CGAL_add_named_parameter(use_compact_clipper_t, use_compact_clipper, use_compact_clipper)
 CGAL_add_named_parameter(output_iterator_t, output_iterator, output_iterator)
+CGAL_add_named_parameter(erase_all_duplicates_t, erase_all_duplicates, erase_all_duplicates)
+CGAL_add_named_parameter(require_same_orientation_t, require_same_orientation, require_same_orientation)
 
 // List of named parameters that we use in the package 'Surface Mesh Simplification'
 CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)

BGL/test/BGL/test_cgal_bgl_named_params.cpp

@@ -85,6 +85,8 @@ void test(const NamedParameters& np)
   assert(get_param(np, CGAL::internal_np::throw_on_self_intersection).v == 43);
   assert(get_param(np, CGAL::internal_np::clip_volume).v == 44);
   assert(get_param(np, CGAL::internal_np::use_compact_clipper).v == 45);
+  assert(get_param(np, CGAL::internal_np::erase_all_duplicates).v == 48);
+  assert(get_param(np, CGAL::internal_np::require_same_orientation).v == 49);
 
     // Named parameters that we use in the package 'Surface Mesh Simplification'
   assert(get_param(np, CGAL::internal_np::get_cost_policy).v == 34);
@@ -162,6 +164,8 @@ void test(const NamedParameters& np)
   check_same_type<43>(get_param(np, CGAL::internal_np::throw_on_self_intersection));
   check_same_type<44>(get_param(np, CGAL::internal_np::clip_volume));
   check_same_type<45>(get_param(np, CGAL::internal_np::use_compact_clipper));
+  check_same_type<48>(get_param(np, CGAL::internal_np::erase_all_duplicates));
+  check_same_type<49>(get_param(np, CGAL::internal_np::require_same_orientation));
 
     // Named parameters that we use in the package 'Surface Mesh Simplification'
   check_same_type<34>(get_param(np, CGAL::internal_np::get_cost_policy));
@@ -241,6 +245,8 @@ int main()
                          .use_compact_clipper(A<45>(45))
                          .apply_per_connected_component(A<46>(46))
                          .output_iterator(A<47>(47))
+                         .erase_all_duplicates(A<48>(48))
+                         .require_same_orientation(A<49>(49))
        );
 
   return EXIT_SUCCESS;

Barycentric_coordinates_2/doc/Barycentric_coordinates_2/Barycentric_coordinates_2.txt

@@ -223,7 +223,7 @@ From the figure above it is easy to see that the \f$O(n^2)\f$ algorithm is as fa
 
 The generic design of the package was developed in 2013 by Dmitry Anisimov and David Bommes with many useful comments by Kai Hormann and Pierre Alliez. The package consists of 6 classes, 2 enumerations, and one namespace. Appropriate iterators are used to provide an efficient access to data and to pass them to one of the generic algorithms for computing coordinates. Once instantiated for a polygon (triangle, segment), the coordinates can be computed multiple times for different query points with respect to all the vertices of the provided polygon (triangle, segment). All the classes are fully templated and have a simple and similar design. In particular, we follow the same naming convention for all functions. Yet, the number of functions can differ from one class to another. 
 
-The implemented algorithms for computing coordinates do not depend on a particular kernel, and all the coordinates can be computed exactly, if an exact kernel is used, apart from mean value coordinates. The latter coordinates involve a square root operation, which results in a slightly worse precision with exact data types due to temporal conversion into a floating point type. The computed coordinates can be stored in an arbitrary container if an appropriate <a href="http://en.cppreference.com/w/cpp/concept/OutputIterator">output iterator</a> is provided. 
+The implemented algorithms for computing coordinates do not depend on a particular kernel, and all the coordinates can be computed exactly, if an exact kernel is used, apart from mean value coordinates. The latter coordinates involve a square root operation, which results in a slightly worse precision with exact data types due to temporal conversion into a floating point type. The computed coordinates can be stored in an arbitrary container if an appropriate <a href="https://en.cppreference.com/w/cpp/named_req/OutputIterator">output iterator</a> is provided.
 
 It is worth noting that the class `CGAL::Barycentric_coordinates::Segment_coordinates_2` is used to compute generalized barycentric coordinates along the polygon's boundary. Hence, one can use the trick for segment coordinates from Section \ref gbc_degeneracies if one is convinced that a point must lie exactly on the polygon's boundary but due to some numerical instabilities it does not.
 

Documentation/doc/Documentation/General.txt

@@ -94,11 +94,37 @@ class RandomAccessIterator {};
 /// See https://en.cppreference.com/w/cpp/named_req/BidirectionalIterator
 class BidirectionalIterator {};
 
+/// \cgalConcept
+/// Concept from the \cpp standard.
+/// See https://en.cppreference.com/w/cpp/named_req/Swappable
+class Swappable {};
+
+/// \cgalConcept
+/// Concept from the \cpp standard.
+/// See https://en.cppreference.com/w/cpp/named_req/Container
+class Container {};
+
+/// \cgalConcept
+/// Concept from the \cpp standard.
+/// See https://en.cppreference.com/w/cpp/named_req/ReversibleContainer
+class ReversibleContainer {};
+
+/// \cgalConcept
+/// Concept from the \cpp standard.
+/// See https://en.cppreference.com/w/cpp/named_req/AssociativeContainer
+class AssociativeContainer {};
+
+/// \cgalConcept
+/// Concept from the \cpp standard.
+/// See https://en.cppreference.com/w/cpp/named_req/SequenceContainer
+class SequenceContainer {};
+
 /// \cgalConcept
 /// This container concept refines
 /// <a href="https://en.cppreference.com/w/cpp/named_req/ReversibleContainer"><tt>ReversibleContainer</tt></a> and its iterator type is a model of
 /// <a href="https://en.cppreference.com/w/cpp/named_req/RandomAccessIterator"><tt>RandomAccessIterator</tt></a>.
 class RandomAccessContainer {};
+
 /// \cgalConcept
 /// This container concepts refines
 /// <a href="https://en.cppreference.com/w/cpp/named_req/SequenceContainer"><tt>SequenceContainer</tt></a> and

Installation/CHANGES.md

@@ -7,7 +7,17 @@ Release 4.14
 Release date: March 2019
 
 ### Polygon Mesh Processing package
--   Added the following new functions to detect and repair mesh degeneracies:
+-   Added the following new functions to detect and repair issues in polygon soups:
+    - `CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup()`, which detects and removes
+       points that are not used in any polygon of the soup.
+    - `CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup()`,
+       which detects and merges points that share the same geometric position.
+    - `CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup()`,
+       which detects and merges polygons that are identical.
+    - `CGAL::Polygon_mesh_processing::repair_polygon_soup()`,
+       which applies a number of repairing steps (a subset of which are the functions above)
+       to clean and repair a polygon soup.
+-   Added the following new functions to detect and repair degeneracies in polygon meshes:
     - `CGAL::Polygon_mesh_processing::degenerate_edges()`
     - `CGAL::Polygon_mesh_processing::degenerate_faces()`
     - `CGAL::Polygon_mesh_processing::is_non_manifold_vertex()`

Mesh_3/include/CGAL/IO/facets_in_complex_3_to_triangle_mesh.h

@@ -54,7 +54,7 @@ void resize(Polygon& p, std::size_t size)
 template <std::size_t N, class INT>
 void resize(CGAL::cpp11::array<INT, N>&, std::size_t CGAL_assertion_code(size))
 {
-  CGAL_assertion(size >= N);
+  CGAL_assertion(size == N);
 }
 
 template<class C3T3, class PointContainer, class FaceContainer>

Mesh_3/include/CGAL/Polyhedral_complex_mesh_domain_3.h

@@ -431,9 +431,9 @@ public:
       needed_vertices_on_patch[i] = (std::min)(nb_of_extra_vertices_per_patch,
                                                needed_vertices_on_patch[i]);
     }
+
     // Then a second path to fill `several_vertices_on_patch`...
-    // The algorithm is adapted from SGI `random_sample_n`:
-    //   https://www.sgi.com/tech/stl/random_sample_n.html
+    // The algorithm is adapted from SGI `random_sample_n`
     BOOST_FOREACH(const Polyhedron& p, this->stored_polyhedra)
     {
       for (typename Polyhedron::Vertex_const_iterator

Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt

@@ -329,15 +329,15 @@ Parameter used in `isotropic_remeshing()` to specify an alternative vertex proje
 \cgalNPBegin{apply_per_connected_component} \anchor PMP_apply_per_connected_component
 Parameter used to indicate whether an algorithm should consider each connected component 
 of a mesh independently.\n
-\b Type : `bool` \n
-\b Default value is `false`
+<b>Type:</b> `bool` \n
+<b>Default:</b> `false`
 \cgalNPEnd
 
 \cgalNPBegin{visitor} \anchor PMP_visitor
 Parameter used to pass a visitor class to a function. Its type and behavior depend on the visited function.
 \n
-\b Type : `A class` \n
-\b Default : Specific to the function visited
+<b>Type:</b> `A class` \n
+<b>Default:</b> Specific to the function visited
 \cgalNPEnd
 
 \cgalNPBegin{throw_on_self_intersection} \anchor PMP_throw_on_self_intersection
@@ -345,23 +345,23 @@ Parameter used in corefinement-related functions to make the functions throw an
 case some faces involved in the intersection of the input are self-intersecting
 and make the operation impossible with the current version of the code.
 \n
-\b Type : `bool` \n
-\b Default value is `false`
+<b>Type:</b> `bool` \n
+<b>Default:</b> `false`
 \cgalNPEnd
 
 \cgalNPBegin{clip_volume} \anchor PMP_clip_volume
 Parameter used in `clip()` functions to clip a volume rather than a surface.
 \n
-\b Type : `bool` \n
-\b Default value is `false`
+<b>Type:</b> `bool` \n
+<b>Default:</b> `false`
 \cgalNPEnd
 
 \cgalNPBegin{use_compact_clipper} \anchor PMP_use_compact_clipper
 Parameter used in `clip()` functions to indicate whether the boundary of the clipper
 should be considered as part of the clipping volume or not.
 \n
-\b Type : `bool` \n
-\b Default value is `true`
+<b>Type:</b> `bool` \n
+<b>Default:</b> `true`
 \cgalNPEnd
 
 \cgalNPBegin{output_iterator} \anchor PMP_output_iterator
@@ -371,6 +371,23 @@ Parameter to pass an output iterator.
 \b Default : `Emptyset_iterator`
 \cgalNPEnd
 
+\cgalNPBegin{erase_all_duplicates} \anchor PMP_erase_all_duplicates
+Parameter used in the function `merge_duplicate_polygons_in_polygon_soup()` to indicate,
+when multiple faces are duplicates, whether all the duplicate faces should be removed
+or if one (arbitrarily chosen) face should be kept.
+\n
+<b>Type:</b> `bool` \n
+<b>Default:</b> `false`
+\cgalNPEnd
+
+\cgalNPBegin{require_same_orientation} \anchor PMP_require_same_orientation
+Parameter used in the function `merge_duplicate_polygons_in_polygon_soup()` to indicate
+if orientation should matter when determining whether two faces are duplicates.
+\n
+<b>Type:</b> `bool` \n
+<b>Default:</b> `false`
+\cgalNPEnd
+
 \cgalNPTableEnd
 
 */

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -36,9 +36,10 @@
 /// \ingroup PkgPolygonMeshProcessingRef
 
 /// \defgroup PMP_repairing_grp Combinatorial Repairing
-/// Functions to orient polygon soups  and to stitch geometrically identical boundaries.
+/// Functions to repair polygon soups and polygon meshes.
 /// \ingroup PkgPolygonMeshProcessingRef
 
+
 /// \defgroup PMP_distance_grp Distance Functions
 /// Functions to compute the distance between meshes, between a mesh and a point set and between a point set and a mesh.
 /// \ingroup PkgPolygonMeshProcessingRef
@@ -128,7 +129,11 @@ and provides a list of the parameters that are used in this package.
 - `CGAL::Polygon_mesh_processing::reverse_face_orientations()`
 
 ## Combinatorial Repairing Functions ##
-- \link PMP_repairing_grp `CGAL::Polygon_mesh_processing::stitch_borders()` \endlink
+- `CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup()`
+- `CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup()`
+- `CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup()`
+- `CGAL::Polygon_mesh_processing::repair_polygon_soup()`
+- `CGAL::Polygon_mesh_processing::stitch_borders()`
 - `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`
 - `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
 - `CGAL::Polygon_mesh_processing::remove_isolated_vertices()`

Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt

@@ -531,7 +531,20 @@ Section \ref PMPOrientation.
 
 ****************************************
 \section PMPRepairing Combinatorial Repairing
+
 *******************
+\subsection PSRepairing Polygon Soup Repairing
+To ensure that a polygon soup can be oriented (see Section \ref PolygonSoups) and transformed
+into a workable polygon mesh, it might be necessary to preprocess the data to remove combinatorial
+and geometrical errors. This package offers the following functions:
+- `CGAL::Polygon_mesh_processing::merge_duplicate_points_in_polygon_soup()`,
+- `CGAL::Polygon_mesh_processing::merge_duplicate_polygons_in_polygon_soup()`,
+- `CGAL::Polygon_mesh_processing::remove_isolated_points_in_polygon_soup()`,
+
+as well as the function `CGAL::Polygon_mesh_processing::repair_polygon_soup()`,
+which bundles the previous functions and an additional handful of repairing techniques
+to obtain an as-clean-as-possible polygon soup.
+
 \subsection Stitching
 
 It happens that a polygon mesh has several edges and vertices that are duplicated.
@@ -575,6 +588,8 @@ is output.
 \endif
 
 \subsection PMPManifoldness Polygon Mesh Manifoldness
+This package offers repairing methods to clean ill-formed polygon soups,
+see Section \ref PMPRepairing.
 
 Non-manifold vertices can be detected using the function `CGAL::Polygon_mesh_processing::is_non_manifold_vertex()`.
 The function `CGAL::Polygon_mesh_processing::duplicate_non_manifold_vertices()` can be used

Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt

@@ -22,4 +22,5 @@
 \example Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp
 \example Polygon_mesh_processing/detect_features_example.cpp
 \example Polygon_mesh_processing/manifoldness_repair_example.cpp
+\example Polygon_mesh_processing/repair_polygon_soup_example.cpp
 */

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -105,6 +105,7 @@ create_single_source_cgal_program( "corefinement_LCC.cpp")
 create_single_source_cgal_program( "hole_filling_example_LCC.cpp" )
 create_single_source_cgal_program( "detect_features_example.cpp" )
 create_single_source_cgal_program( "manifoldness_repair_example.cpp" )
+create_single_source_cgal_program( "repair_polygon_soup_example.cpp" )
 
 if(OpenMesh_FOUND)
 create_single_source_cgal_program( "compute_normals_example_OM.cpp" )

Polygon_mesh_processing/examples/Polygon_mesh_processing/repair_polygon_soup_example.cpp

@@ -0,0 +1,82 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+
+#include <iostream>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
+typedef K::Point_3                                              Point_3;
+
+typedef std::vector<std::size_t>                                Polygon;
+typedef CGAL::Surface_mesh<Point_3>                             Mesh;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(int, char**)
+{
+  // First, construct a polygon soup with some problems
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  points.push_back(Point_3(0,0,0));
+  points.push_back(Point_3(1,0,0));
+  points.push_back(Point_3(0,1,0));
+  points.push_back(Point_3(-1,0,0));
+  points.push_back(Point_3(0,-1,0));
+  points.push_back(Point_3(0,1,0)); // duplicate point
+  points.push_back(Point_3(0,-2,0)); // unused point
+
+  Polygon p;
+  p.push_back(0); p.push_back(1); p.push_back(2);
+  polygons.push_back(p);
+
+  // degenerate face
+  p.clear();
+  p.push_back(0); p.push_back(0); p.push_back(0);
+  polygons.push_back(p);
+
+  p.clear();
+  p.push_back(0); p.push_back(1); p.push_back(4);
+  polygons.push_back(p);
+
+  // duplicate face with different orientation
+  p.clear();
+  p.push_back(0); p.push_back(4); p.push_back(1);
+  polygons.push_back(p);
+
+  p.clear();
+  p.push_back(0); p.push_back(3); p.push_back(5);
+  polygons.push_back(p);
+
+  // degenerate face
+  p.clear();
+  p.push_back(0); p.push_back(3); p.push_back(0);
+  polygons.push_back(p);
+
+  p.clear();
+  p.push_back(0); p.push_back(3); p.push_back(4);
+  polygons.push_back(p);
+
+  // pinched and degenerate face
+  p.clear();
+  p.push_back(0); p.push_back(1); p.push_back(2); p.push_back(3);
+  p.push_back(4); p.push_back(3); p.push_back(2); p.push_back(1);
+  polygons.push_back(p);
+
+  PMP::repair_polygon_soup(points, polygons);
+  PMP::orient_polygon_soup(points, polygons);
+
+  Mesh mesh;
+  PMP::polygon_soup_to_polygon_mesh(points, polygons, mesh);
+
+  std::cout << "Mesh has " << num_vertices(mesh) << " vertices and " << num_faces(mesh) << " faces" << std::endl;
+
+  assert(num_vertices(mesh) == 5);
+  assert(num_faces(mesh) == 4);
+
+  return 0;
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/connected_components.h

@@ -277,7 +277,7 @@ void keep_connected_components(PolygonMesh& pmesh
  *
  * \param pmesh the polygon mesh
  * \param nb_components_to_keep the number of components to be kept
- * \param np optional \ref pmp_namedparameters "Named Parameters" described below
+ * \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
  *
  * \cgalNamedParamsBegin
  *    \cgalParamBegin{edge_is_constrained_map} a property map containing the constrained-or-not status of each edge of `pmesh` \cgalParamEnd
@@ -360,7 +360,7 @@ std::size_t keep_largest_connected_components(PolygonMesh& pmesh,
  *
  * \param pmesh the polygon mesh
  * \param threshold_components_to_keep the number of faces a component must have so that it is kept
- * \param np optional \ref pmp_namedparameters "Named Parameters" described below
+ * \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
  *
  * \cgalNamedParamsBegin
  *    \cgalParamBegin{edge_is_constrained_map} a property map containing the constrained-or-not status of each edge of `pmesh` \cgalParamEnd
@@ -664,7 +664,7 @@ void remove_connected_components(PolygonMesh& pmesh
 *
 * \param components_to_remove a face range, including one face or more on each component to be removed
 * \param pmesh the polygon mesh
-* \param np optional \ref pmp_namedparameters "Named Parameters" described below
+* \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{edge_is_constrained_map} a property map containing the constrained-or-not status of each edge of `pmesh` \cgalParamEnd
@@ -722,7 +722,7 @@ void remove_connected_components(PolygonMesh& pmesh
 *
 * \param pmesh the polygon mesh
 * \param components_to_keep a face range, including one face or more on each component to be kept
-* \param np optional \ref pmp_namedparameters "Named Parameters" described below
+* \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
 *
 * \cgalNamedParamsBegin
 *    \cgalParamBegin{edge_is_constrained_map} a property map containing the constrained-or-not status of each edge of `pmesh` \cgalParamEnd

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/detect_features.h

@@ -260,7 +260,7 @@ template<typename GT,
  * \param pmesh the polygon mesh
  * \param angle_in_deg the dihedral angle bound
  * \param edge_is_feature_map the property map that will contain the sharp-or-not status of each edge of `pmesh`
- * \param np optional \ref pmp_namedparameters "Named Parameters" described below
+ * \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
  *
  * \cgalNamedParamsBegin
  *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd
@@ -400,7 +400,7 @@ namespace internal
  * \param angle_in_deg the dihedral angle bound
  * \param edge_is_feature_map the property map that will contain the sharp-or-not status of each edge of `pmesh`
  * \param patch_id_map the property map that will contain the surface patch ids for the faces of `pmesh`.
- * \param np optional \ref pmp_namedparameters "Named Parameters" described below
+ * \param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
  *
  * \cgalNamedParamsBegin
  *    \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/orient_polygon_soup.h

@@ -463,13 +463,13 @@ struct Polygon_soup_orienter
  * The algorithm is described in \cgalCite{gueziec2001cutting}.
  *
  * @tparam PointRange a model of the concepts `RandomAccessContainer`
- * and `BackInsertionSequence` whose value type is the point type
+ * and `BackInsertionSequence` whose value type is the point type.
  * @tparam PolygonRange a model of the concept `RandomAccessContainer`
  * whose value_type is a model of the concept `RandomAccessContainer`
  * whose value_type is `std::size_t`.
  *
- * @param points points of the soup of polygons. Some points might be pushed back to resolve
- *               non-manifold or non-orientability issues.
+ * @param points points of the soup of polygons. Some additional points might be pushed back to resolve
+ *               non-manifoldness or non-orientability issues.
  * @param polygons each element in the vector describes a polygon using the index of the points in `points`.
  *                 If needed the order of the indices of a polygon might be reversed.
  * @return `true`  if the orientation operation succeded.

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h

@@ -204,7 +204,8 @@ public:
   * @param polygons each element in the vector describes a polygon using the index of the points in `points`
   * @param out the polygon mesh to be built
   *
-  * @pre `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(polygons)`
+  * @pre \link CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()
+  *            CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(polygons) \endlink
   *
   * \sa `CGAL::Polygon_mesh_processing::orient_polygon_soup()`
   * \sa `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h

@@ -309,7 +309,7 @@ void isotropic_remeshing(
 * @param edges the range of edges to be split if they are longer than given threshold
 * @param max_length the edge length above which an edge from `edges` is split
 *        into to sub-edges
-* @param np optional \ref pmp_namedparameters "Named Parameters" described below
+* @param np optional \ref pmp_namedparameters "Named Parameters", amongst those described below
 
 * \cgalNamedParamsBegin
 *  \cgalParamBegin{vertex_point_map} the property map with the points associated

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair.h

@@ -850,12 +850,12 @@ bool remove_degenerate_edges(const EdgeRange& edge_range,
 //    \cgalParamEnd
 //    \cgalParamBegin{geom_traits} a geometric traits class instance.
 //       The traits class must provide the nested type `Point_3`,
-//       and the nested functors :
+//       and the nested functors:
 //         - `Compare_distance_3` to compute the distance between 2 points
 //         - `Collinear_3` to check whether 3 points are collinear
 //         - `Less_xyz_3` to compare lexicographically two points
-//         - `Equal_3` to check whether 2 points are identical
-//       For each functor `Foo`, a function `Foo foo_object()`
+///        - `Equal_3` to check whether 2 points are identical.
+///       For each functor Foo, a function `Foo foo_object()` must be provided.
 //   \cgalParamEnd
 // \cgalNamedParamsEnd
 //

Polygon_mesh_processing/include/CGAL/polygon_mesh_processing.h

@@ -50,6 +50,7 @@
 #include <CGAL/Polygon_mesh_processing/intersection.h>
 #include <CGAL/Polygon_mesh_processing/transform.h>
 #include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
 #include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
 
 // the named parameter header being not documented the doc is put here for now

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -102,6 +102,7 @@ endif()
   create_single_source_cgal_program("test_pmp_transform.cpp")
   create_single_source_cgal_program("extrude_test.cpp")
   create_single_source_cgal_program("test_merging_border_vertices.cpp")
+  create_single_source_cgal_program("test_repair_polygon_soup.cpp")
   create_single_source_cgal_program("test_shape_predicates.cpp")
   create_single_source_cgal_program("test_pmp_collision_detection.cpp")
 

Polygon_mesh_processing/test/Polygon_mesh_processing/test_repair_polygon_soup.cpp

@@ -0,0 +1,551 @@
+#define CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE_PP
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
+
+#include <CGAL/Surface_mesh.h>
+
+#include <algorithm>
+#include <deque>
+#include <iostream>
+#include <fstream>
+#include <vector>
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+namespace params = PMP::parameters;
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
+typedef K::Point_3                                              Point_3;
+
+typedef CGAL::Surface_mesh<Point_3>                             Mesh;
+typedef std::vector<std::size_t>                                Polygon;
+
+void test_polygon_canonicalization(const bool verbose = false)
+{
+  std::cout << "test polygon canonicalization... " << std::endl;
+
+  std::vector<Point_3> points;
+  points.push_back(Point_3(0,0,0)); // #0
+  points.push_back(Point_3(1,0,0)); // #1
+  points.push_back(Point_3(0,1,0)); // #2
+  points.push_back(Point_3(1,1,0)); // #3
+  points.push_back(Point_3(1,1,2)); // #4
+  points.push_back(Point_3(1,1,-2)); // #5
+
+  // empty
+  Polygon polygon;
+  Polygon canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon.empty());
+
+  // 1 point
+  polygon.push_back(5);
+
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon == polygon);
+
+  // 2 points
+  polygon.clear();
+  polygon.push_back(4); polygon.push_back(3);
+
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon[0] == 3 && canonical_polygon[1] == 4);
+  std::swap(polygon[0], polygon[1]);
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon[0] == 3 && canonical_polygon[1] == 4);
+
+  // 3 points
+  polygon.clear();
+  polygon.push_back(4); polygon.push_back(1); polygon.push_back(3);
+
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon[0] == 1 && canonical_polygon[1] == 3 && canonical_polygon[2] == 4);
+  std::swap(polygon[0], polygon[2]);
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon[0] == 1 && canonical_polygon[1] == 3 && canonical_polygon[2] == 4);
+
+  // Generic case
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(5); polygon.push_back(4); polygon.push_back(1);
+
+  // Whether reversed or with cyclic permutations, it should always yield the same canonical polygon
+  canonical_polygon = PMP::internal::construct_canonical_polygon(points, polygon, K());
+  assert(canonical_polygon.size() == 5);
+  assert(canonical_polygon[0] == 0);
+
+  // all cyclic permutations
+  for(std::size_t i=0, end=polygon.size(); i<end; ++i)
+  {
+    Polygon cpol = PMP::internal::construct_canonical_polygon(points, polygon, K());
+    if(verbose)
+    {
+      std::cout << "Input polygon:";
+      PMP::internal::print_polygon(std::cout, polygon);
+      std::cout << "Canonical polygon:";
+      PMP::internal::print_polygon(std::cout, canonical_polygon);
+    }
+    assert(cpol == canonical_polygon);
+
+    std::rotate(polygon.begin(), polygon.begin()+1, polygon.end());
+  }
+
+  // reverse and all cyclic permutations
+  std::reverse(polygon.begin(), polygon.end());
+  for(std::size_t i=0, end=polygon.size(); i<end; ++i)
+  {
+    Polygon cpol = PMP::internal::construct_canonical_polygon(points, polygon, K());
+    if(verbose)
+    {
+      std::cout << "Input polygon:";
+      PMP::internal::print_polygon(std::cout, polygon);
+      std::cout << "Canonical polygon:";
+      PMP::internal::print_polygon(std::cout, canonical_polygon);
+    }
+    assert(cpol == canonical_polygon);
+
+    std::rotate(polygon.begin(), polygon.begin()+1, polygon.end());
+  }
+}
+
+void test_merge_duplicate_points(const bool /*verbose*/ = false)
+{
+  std::cout << "test merge duplicate points... " << std::endl;
+
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  // empty
+  std::size_t res = PMP::merge_duplicate_points_in_polygon_soup(points, polygons);
+  assert(res == 0 && points.empty() && polygons.empty());
+
+  points.push_back(Point_3(0,0,0)); // #0
+  points.push_back(Point_3(1,2,0)); // #1
+  points.push_back(Point_3(1,1,0)); // #2
+  points.push_back(Point_3(1,1,0)); // #3 // identical to #2
+  points.push_back(Point_3(0,1,0)); // #4
+  points.push_back(Point_3(1,1,0)); // #5 // identical to #2
+  points.push_back(Point_3(0,0,0)); // #6 // idental to #0
+
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(1); polygon.push_back(2);
+  polygons.push_back(polygon);
+
+  polygon.clear();
+  polygon.push_back(6); polygon.push_back(3); polygon.push_back(1); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  polygon.clear();
+  polygon.push_back(5);
+  polygons.push_back(polygon);
+
+  res = PMP::merge_duplicate_points_in_polygon_soup(points, polygons, params::geom_traits(K()));
+  assert(res == 3 && points.size() == 4 && polygons.size() == 3);
+
+  assert(polygons[0][0] == 0 && polygons[0][1] == 1 && polygons[0][2] == 2);
+  assert(polygons[1][0] == 0 && polygons[1][1] == 2 && polygons[1][2] == 1 && polygons[1][3] == 0);
+
+  for(std::size_t i=0, psn=polygons.size(); i<psn; ++i)
+  {
+    const Polygon& polygon = polygons[i];
+    for(std::size_t j=0, pn=polygon.size(); j<pn; ++j)
+    {
+      assert(polygon[j] < points.size());
+    }
+  }
+}
+
+void test_merge_duplicate_polygons(const bool /*verbose*/ = false)
+{
+  std::cout << "test duplicate polygons merging..." << std::endl;
+
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  points.push_back(Point_3(0,0,0)); // #0
+  points.push_back(Point_3(1,0,0)); // #1
+  points.push_back(Point_3(0,1,0)); // #2
+  points.push_back(Point_3(1,1,0)); // #3
+  points.push_back(Point_3(1,1,2)); // #4
+
+  // -------------------------------------------------------
+  // empty
+  std::vector<std::vector<std::size_t> > all_duplicate_polygons;
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), true /*only equal if same orientation*/);
+  assert(polygons.empty() && all_duplicate_polygons.empty());
+
+  std::size_t res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons, params::geom_traits(K()));
+  assert(res == 0 && polygons.empty());
+
+  // -------------------------------------------------------
+  // 1 polygon
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(1); polygon.push_back(2);
+  polygons.push_back(polygon);
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), false /*equal regardless of orientation*/);
+  assert(all_duplicate_polygons.empty());
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), true /*only equal if same orientation*/);
+  assert(all_duplicate_polygons.empty());
+
+  res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons, params::geom_traits(K()));
+  assert(res == 0 && polygons.size() == 1);
+
+  // -------------------------------------------------------
+  // 2 different polygons
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(1); polygon.push_back(3); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), false /*equal regardless of orientation*/);
+  assert(all_duplicate_polygons.empty());
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), true /*only equal if same orientation*/);
+  assert(all_duplicate_polygons.empty());
+
+  res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons);
+  assert(res == 0 && polygons.size() == 2);
+
+  // -------------------------------------------------------
+  // Multiple duplicates
+    // duplicate, same orientations
+  polygon.clear();
+  polygon.push_back(2); polygon.push_back(0); polygon.push_back(1);
+  polygons.push_back(polygon);
+
+    // duplicate, different orientations
+  polygon.clear();
+  polygon.push_back(2); polygon.push_back(1); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+    // duplicate, different orientations
+  polygon.clear();
+  polygon.push_back(4); polygon.push_back(3); polygon.push_back(1); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+    // same vertices, not a duplicate
+  polygon.clear();
+  polygon.push_back(3); polygon.push_back(4); polygon.push_back(1); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), true /*only equal if same orientation*/);
+  assert(all_duplicate_polygons.size() == 1); // one duplication
+  assert(all_duplicate_polygons[0].size() == 2); // two polygons are equal
+  all_duplicate_polygons.clear();
+
+  PMP::internal::collect_duplicate_polygons(points, polygons,
+                                            std::back_inserter(all_duplicate_polygons),
+                                            K(), false /*equal regardless of orientation*/);
+  assert(all_duplicate_polygons.size() == 2);
+
+  // not sure which duplicate is output first
+  if(all_duplicate_polygons[0][0] == 0)
+    assert(all_duplicate_polygons[0].size() == 3 && all_duplicate_polygons[1].size() == 2);
+  else
+    assert(all_duplicate_polygons[0].size() == 2 && all_duplicate_polygons[1].size() == 3);
+
+  // Keep one for each duplicate
+  std::vector<Polygon> polygons_copy(polygons);
+  res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons_copy,
+                                                      params::all_default());
+  assert(res == 3 && polygons_copy.size() == 3);
+
+  // Remove all duplicates
+  polygons_copy = polygons;
+  res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons_copy,
+                                                      params::erase_all_duplicates(true)
+                                                             .require_same_orientation(false));
+  assert(res == 5 && polygons_copy.size() == 1);
+
+  // Remove all duplicates but different orientations are different polygons
+  res = PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons,
+                                                      params::erase_all_duplicates(true)
+                                                             .require_same_orientation(true));
+  assert(res == 2 && polygons.size() == 4);
+}
+
+void test_simplify_polygons(const bool /*verbose*/ = false)
+{
+  std::cout << "test simplify_polygons... " << std::endl;
+
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  points.push_back(Point_3(0,0,0)); // #0
+  points.push_back(Point_3(1,2,0)); // #1
+  points.push_back(Point_3(1,0,0)); // #2
+  points.push_back(Point_3(1,3,0)); // #3
+  points.push_back(Point_3(0,1,0)); // #4
+  points.push_back(Point_3(1,1,0)); // #5
+  points.push_back(Point_3(0,0,0)); // #6
+
+  // ------
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  std::size_t res = PMP::internal::simplify_polygons_in_polygon_soup<K>(points, polygons);
+  assert(res == 0 && polygons.back().size() == 3);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 1);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(0); polygon.push_back(0); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 1);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(3); polygon.push_back(3); polygon.push_back(3); polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 2);
+
+  // ------
+  // Now with the same geometric positions, but different combinatorial information
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(1); polygon.push_back(6);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 3);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(6); polygon.push_back(0); polygon.push_back(5);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 2);
+
+  // ------
+  polygon.clear();
+  polygon.push_back(0); polygon.push_back(6); polygon.push_back(5); polygon.push_back(3); polygon.push_back(3);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::simplify_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.back().size() == 3);
+}
+
+void test_remove_invalid_polygons(const bool /*verbose*/ = false)
+{
+  std::cout << "test remove_invalid_polygons... " << std::endl;
+
+  // points are not actually needed since only the size of the polygons is considered
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  std::size_t res = PMP::internal::remove_invalid_polygons_in_polygon_soup(points, polygons);
+  assert(res == 0 && polygons.size() == 0);
+
+  // non-trivial polygon
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::remove_invalid_polygons_in_polygon_soup(points, polygons);
+  assert(res == 0 && polygons.size() == 1);
+
+  // another non-trivial polygon
+  polygon.clear();
+  polygon.push_back(1); polygon.push_back(3); polygon.push_back(3); polygon.push_back(7);
+  polygons.push_back(polygon);
+
+  // empty polygon
+  polygon.clear();
+  polygons.push_back(polygon);
+
+  // 1-vertex polygon
+  polygon.push_back(0);
+  polygons.push_back(polygon);
+
+  // 2-vertex polygon
+  polygon.push_back(1);
+  polygons.push_back(polygon);
+
+  // another non-trivial polygon
+  polygon.clear();
+  polygon.push_back(8); polygon.push_back(9); polygon.push_back(7); polygon.push_back(6);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::remove_invalid_polygons_in_polygon_soup(points, polygons);
+  assert(res == 3 && polygons.size() == 3);
+}
+
+template <typename PointRange, typename PolygonRange>
+std::size_t test_remove_isolated_points_data_set(PointRange& points,
+                                                 PolygonRange& polygons,
+                                                 const bool verbose = false)
+{
+  if(verbose)
+  {
+    std::cout << "Input:" << std::endl;
+    std::cout << points.size() << " points" << std::endl;
+    std::cout << polygons.size() << " polygons" << std::endl;
+  }
+
+  std::size_t rm_nv = PMP::remove_isolated_points_in_polygon_soup(points, polygons);
+
+  if(verbose)
+  {
+    std::cout << "Removed " << rm_nv << " points" << std::endl;
+    std::cout << points.size() << " points" << std::endl;
+    std::cout << polygons.size() << " polygons" << std::endl;
+
+    std::cout << "Polygons:" << std::endl;
+    for(std::size_t i=0; i<polygons.size(); ++i)
+    {
+      for(std::size_t j=0; j<polygons[i].size(); ++j)
+        std::cout << polygons[i][j] << " ";
+      std::cout << std::endl;
+    }
+  }
+
+  return rm_nv;
+}
+
+void test_remove_isolated_points(const bool verbose = false)
+{
+  std::cout << "test remove_isolated_points... " << std::endl;
+
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  // everything empty
+  std::size_t res = test_remove_isolated_points_data_set(points, polygons, verbose);
+  assert(res == 0 && points.empty() && polygons.empty());
+
+  points.push_back(Point_3(0,0,0));
+  points.push_back(Point_3(1,2,0));
+  points.push_back(Point_3(1,0,0));
+  points.push_back(Point_3(1,3,0));
+  points.push_back(Point_3(0,1,0));
+  points.push_back(Point_3(1,1,0));
+
+  // no polygons (all points are unused)
+  std::deque<Point_3> points_copy(points.begin(), points.end());
+  res = test_remove_isolated_points_data_set(points_copy, polygons, verbose);
+  assert(res == 6 && points_copy.empty() && polygons.empty());
+
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  polygon.clear();
+  polygon.push_back(4); polygon.push_back(0); polygon.push_back(2); polygon.push_back(5);
+  polygons.push_back(polygon);
+
+  // generic test
+  res = test_remove_isolated_points_data_set(points, polygons, verbose);
+  assert(res == 2 && points.size() == 4 && polygons.size() == 2);
+}
+
+void test_slit_pinched_polygons(const bool /*verbose*/ = false)
+{
+  std::cout << "test split_pinched_polygons... " << std::endl;
+
+  std::vector<Point_3> points;
+  std::vector<Polygon> polygons;
+
+  // everything empty
+  std::size_t res = PMP::internal::split_pinched_polygons_in_polygon_soup<K>(points, polygons);
+  assert(res == 0 && points.empty() && polygons.empty());
+
+  points.push_back(Point_3(0,0,0)); // #0
+  points.push_back(Point_3(1,0,0)); // #1
+  points.push_back(Point_3(0,1,0)); // #2
+  points.push_back(Point_3(1,1,0)); // #3
+  points.push_back(Point_3(1,3,0)); // #4
+  points.push_back(Point_3(1,1,0)); // #5
+
+  // no pinch
+  Polygon polygon;
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::split_pinched_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 0 && polygons.size() == 1);
+
+  // pinch via same ID (2)
+  polygon.clear();
+  polygons.clear();
+  polygon.push_back(0); polygon.push_back(2); polygon.push_back(4); polygon.push_back(3); polygon.push_back(2);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::split_pinched_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.size() == 2);
+  assert(polygons[0].size() == 3 && polygons[1].size() == 2);
+
+  // pinch via same point (5 & 3)
+  polygon.clear();
+  polygons.clear();
+  polygon.push_back(5); polygon.push_back(1); polygon.push_back(0);
+  polygon.push_back(3); polygon.push_back(4); polygon.push_back(2);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::split_pinched_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.size() == 2);
+  assert(polygons[0].size() == 3 && polygons[1].size() == 3);
+
+  // pinch on last point
+  polygon.clear();
+  polygons.clear();
+  polygon.push_back(1); polygon.push_back(5); polygon.push_back(0);
+  polygon.push_back(2); polygon.push_back(4); polygon.push_back(3);
+  polygons.push_back(polygon);
+
+  res = PMP::internal::split_pinched_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 1 && polygons.size() == 2);
+  assert(polygons[0].size() == 4 && polygons[1].size() == 2);
+
+  // multiple pinches
+  polygon.clear();
+  polygons.clear();
+  polygon.push_back(5); polygon.push_back(1); polygon.push_back(3); // pinch 5&3, pinch 3...
+  polygon.push_back(2); polygon.push_back(4); polygon.push_back(0);
+  polygon.push_back(1); polygon.push_back(3); polygon.push_back(1); // ... and 3, pinch 1...
+  polygon.push_back(2); polygon.push_back(4); polygon.push_back(1); // ... and 1
+  polygons.push_back(polygon);
+
+  res = PMP::internal::split_pinched_polygons_in_polygon_soup(points, polygons, K());
+  assert(res == 3 && polygons.size() == 4);
+  assert(polygons[0].size() == 2); // 5 1
+  assert(polygons[1].size() == 5); // 3 2 4 5 1
+  assert(polygons[2].size() == 3); // 1 2 4
+  assert(polygons[3].size() == 2); // 1 3
+}
+
+int main()
+{
+  test_polygon_canonicalization(true);
+  test_merge_duplicate_points(false);
+  test_merge_duplicate_polygons(false);
+  test_simplify_polygons(false);
+  test_remove_invalid_polygons(false);
+  test_remove_isolated_points(false);
+  test_slit_pinched_polygons(false);
+
+  return EXIT_SUCCESS;
+}

Polyhedron_IO/include/CGAL/IO/OFF_reader.h

@@ -21,6 +21,7 @@
 #define CGAL_IO_OFF_READER_H
 
 #include <CGAL/IO/File_scanner_OFF.h>
+#include <CGAL/IO/reader_helpers.h>
 
 #include <vector>
 #include <iostream>
@@ -28,33 +29,7 @@
 #include <CGAL/assertions.h>
 #include <CGAL/use.h>
 
-namespace CGAL{
-
-  namespace read_OFF_internal{
-    template <class Point_3>
-    void fill_point(double x, double y, double z, Point_3& pt)
-    {
-      pt = Point_3(x, y, z);
-    }
-
-    inline void fill_point(double x, double y, double z, CGAL::cpp11::array<double,3>& p)
-    {
-      p = CGAL::make_array(x,y,z);
-    }
-
-    template <class Polygon_3>
-    void resize(Polygon_3& p, std::size_t size)
-    {
-      p.resize(size);
-    }
-
-    template <std::size_t N, class INT>
-    void resize(CGAL::cpp11::array<INT, N>&, std::size_t size)
-    {
-      CGAL_USE(size);
-      CGAL_assertion( size>=N );
-    }
-  }
+namespace CGAL {
 
   template <class Point_3, class Polygon_3>
   bool
@@ -71,7 +46,7 @@ namespace CGAL{
       double x, y, z, w;
       scanner.scan_vertex( x, y, z, w);
       CGAL_assertion(w!=0);
-      read_OFF_internal::fill_point( x/w, y/w, z/w, points[i] );
+      IO::internal::fill_point( x/w, y/w, z/w, points[i] );
       scanner.skip_to_next_vertex( i);
     }
     if(!in)
@@ -81,7 +56,7 @@ namespace CGAL{
       std::size_t no;
 
       scanner.scan_facet( no, i);
-      read_OFF_internal::resize(polygons[i], no);
+      IO::internal::resize(polygons[i], no);
       for(std::size_t j = 0; j < no; ++j) {
         std::size_t id;
         scanner.scan_facet_vertex_index(id, i);
@@ -114,7 +89,7 @@ namespace CGAL{
         double x, y, z, w;
         scanner.scan_vertex( x, y, z, w);
         CGAL_assertion(w!=0);
-        read_OFF_internal::fill_point( x/w, y/w, z/w, points[i] );
+        IO::internal::fill_point( x/w, y/w, z/w, points[i] );
         if(scanner.has_colors())
         {
             unsigned char r=0, g=0, b=0;
@@ -131,7 +106,7 @@ namespace CGAL{
       std::size_t no;
       scanner.scan_facet( no, i);
 
-      read_OFF_internal::resize(polygons[i], no);
+      IO::internal::resize(polygons[i], no);
       for(std::size_t j = 0; j < no; ++j) {
         std::size_t id;
         scanner.scan_facet_vertex_index(id, i);

Polyhedron_IO/include/CGAL/IO/STL_reader.h

@@ -21,8 +21,8 @@
 #ifndef CGAL_IO_STL_READER_H
 #define CGAL_IO_STL_READER_H
 
-#include <CGAL/array.h>
 #include <CGAL/IO/io.h>
+#include <CGAL/IO/reader_helpers.h>
 
 #include <boost/cstdint.hpp> 
 
@@ -34,11 +34,12 @@
 
 namespace CGAL {
 
+template <class Point, class Triangle>
 bool read_ASCII_facet(std::istream& input,
-                      std::vector<cpp11::array<double,3> >& points,
-                      std::vector<cpp11::array<int,3> >& facets,
+                      std::vector<Point>& points,
+                      std::vector<Triangle>& facets,
                       int& index,
-                      std::map<cpp11::array<double, 3>, int >& index_map,
+                      std::map<Point, int >& index_map,
                       bool verbose = false)
 {
   // Here, we have already read the word 'facet' and are looking to read till 'endfacet'
@@ -48,8 +49,10 @@ bool read_ASCII_facet(std::istream& input,
               endfacet("endfacet");
 
   int count = 0;
-  cpp11::array<double, 3> p;
-  cpp11::array<int, 3> ijk;
+  double x,y,z;
+  Point p;
+  Triangle ijk;
+  IO::internal::resize(ijk, 3);
 
   while(input >> s)
   {
@@ -76,7 +79,7 @@ bool read_ASCII_facet(std::istream& input,
         return false;
       }
 
-      if(!(input >> iformat(p[0]) >> iformat(p[1]) >> iformat(p[2])))
+      if(!(input >> iformat(x) >> iformat(y) >> iformat(z)))
       {
         if(verbose)
           std::cerr << "Error while reading point coordinates (premature end of file)" << std::endl;
@@ -85,8 +88,8 @@ bool read_ASCII_facet(std::istream& input,
       }
       else
       {
-        std::map<cpp11::array<double, 3>, int>::iterator iti=
-            index_map.insert(std::make_pair(p, -1)).first;
+        IO::internal::fill_point(x, y, z, p);
+        typename std::map<Point, int>::iterator iti = index_map.insert(std::make_pair(p, -1)).first;
 
         if(iti->second == -1)
         {
@@ -110,9 +113,10 @@ bool read_ASCII_facet(std::istream& input,
   return false;
 }
 
+template <class Point, class Triangle>
 bool parse_ASCII_STL(std::istream& input,
-                     std::vector<cpp11::array<double,3> >& points,
-                     std::vector<cpp11::array<int,3> >& facets,
+                     std::vector<Point>& points,
+                     std::vector<Triangle>& facets,
                      bool verbose = false)
 {
   if(verbose)
@@ -124,11 +128,9 @@ bool parse_ASCII_STL(std::istream& input,
   // Here, we have already read the word 'solid'
 
   int index = 0;
-  std::map<cpp11::array<double, 3>, int> index_map;
+  std::map<Point, int> index_map;
 
-  std::string s;
-  std::string facet("facet"),
-              endsolid("endsolid");
+  std::string s, facet("facet"), endsolid("endsolid");
 
   while(input >> s)
   {
@@ -149,9 +151,10 @@ bool parse_ASCII_STL(std::istream& input,
   return false;
 }
 
+template <typename Point, typename Triangle>
 bool parse_binary_STL(std::istream& input,
-                      std::vector<cpp11::array<double,3> >& points,
-                      std::vector<cpp11::array<int,3> >& facets,
+                      std::vector<Point>& points,
+                      std::vector<Triangle>& facets,
                       bool verbose = false)
 {
   if(verbose)
@@ -190,7 +193,7 @@ bool parse_binary_STL(std::istream& input,
     return true; // empty file
 
   int index = 0;
-  std::map<cpp11::array<double, 3>, int> index_map;
+  std::map<Point, int> index_map;
 
   boost::uint32_t N32;
   if(!(input.read(reinterpret_cast<char*>(&N32), sizeof(N32))))
@@ -218,7 +221,8 @@ bool parse_binary_STL(std::istream& input,
       return false;
     }
 
-    cpp11::array<int, 3> ijk;
+    Triangle ijk;
+    IO::internal::resize(ijk, 3);
 
     for(int j=0; j<3; ++j)
     {
@@ -233,11 +237,10 @@ bool parse_binary_STL(std::istream& input,
         return false;
       }
 
-      cpp11::array<double, 3> p;
-      p[0] = x; p[1] = y; p[2] = z;
+      Point p;
+      IO::internal::fill_point(x, y, z, p);
 
-      std::map<cpp11::array<double, 3>, int>::iterator iti =
-        index_map.insert(std::make_pair(p, -1)).first;
+      typename std::map<Point, int>::iterator iti = index_map.insert(std::make_pair(p, -1)).first;
 
       if(iti->second == -1)
       {
@@ -268,9 +271,29 @@ bool parse_binary_STL(std::istream& input,
   return true;
 }
 
+//
+// Read a file with `.stl` format.
+//
+// \tparam Point must be a model of the concept `RandomAccessContainer` or a %CGAL point type
+// \tparam Triangle must be a model of the concept `RandomAccessContainer`
+//
+// \param input the input stream
+// \param points a container that will contain the points used in the .stl file
+// \param polygons a container that will contain the triangles used in the .stl file
+// \param verbose whether to enable or not a sanity log
+//
+// \returns `true` if the reading process went well, `false` otherwise
+//
+// \warning `points` and `facets` are not cleared: new points and triangles are added to the back
+//          of the containers.
+//
+// Although the STL file format uses triangles, it is convenient to be able to use vectors
+// and other models of the `SequenceContainer` (instead of arrays) for the face type,
+// to avoid having to convert the to apply polygon soup reparation algorithms.
+template <class Point, class Triangle>
 bool read_STL(std::istream& input,
-              std::vector<cpp11::array<double,3> >& points,
-              std::vector<cpp11::array<int,3> >& facets,
+              std::vector<Point>& points,
+              std::vector<Triangle>& facets,
               bool verbose = false)
 {
   int pos = 0;

Polyhedron_IO/include/CGAL/IO/reader_helpers.h

@@ -0,0 +1,90 @@
+// Copyright (c) 2015 GeometryFactory
+//
+// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation; either version 3 of the License,
+// or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0+
+//
+// Author(s)     : Mael Rouxel-Labbé
+
+#ifndef CGAL_IO_READER_HELPERS_H
+#define CGAL_IO_READER_HELPERS_H
+
+#include <CGAL/array.h>
+#include <CGAL/assertions.h>
+#include <CGAL/Has_member.h>
+
+#include <boost/mpl/logical.hpp>
+#include <boost/utility/enable_if.hpp>
+
+namespace CGAL{
+
+namespace IO {
+
+namespace internal {
+
+CGAL_GENERATE_MEMBER_DETECTOR(size);
+CGAL_GENERATE_MEMBER_DETECTOR(resize);
+
+// Typical container
+template <class Container>
+void resize(Container& c, std::size_t size,
+            typename boost::enable_if_c<has_resize<Container>::value>::type* = NULL)
+{
+  c.resize(size);
+}
+
+// Container without a resize() function, but with a size() function (e.g. an array)
+template <class Container>
+void resize(Container& CGAL_assertion_code(array), std::size_t CGAL_assertion_code(size),
+            typename boost::enable_if<
+              boost::mpl::and_<
+                boost::mpl::not_<has_resize<Container> >,
+                                 has_size<Container> > >::type* = NULL)
+{
+  CGAL_assertion(array.size() == size);
+}
+
+// A class with neither resize() nor size(), can't enforce size (it better be correct!)
+template <class Container>
+void resize(Container&, std::size_t,
+            typename boost::disable_if<
+              boost::mpl::or_<has_resize<Container>,
+                              has_size<Container> > >::type* = NULL)
+{
+}
+
+// Ideally this should be a std::is_constructible(double, double, double) but boost::is_constructible
+// is not safe to use without CXX11
+template <typename Kernel>
+void fill_point(const double x, const double y, const double z, CGAL::Point_3<Kernel>& pt)
+{
+  pt = CGAL::Point_3<Kernel>(x, y, z);
+}
+
+template <typename Point_3>
+void fill_point(const double x, const double y, const double z, Point_3& pt)
+{
+  // just in case something weirder than arrays or CGAL points are used as points...
+  resize(pt, 3);
+
+  pt[0] = x; pt[1] = y; pt[2] = z;
+}
+
+} // end namespace internal
+
+} // end namespace IO
+
+} // namespace CGAL
+
+#endif // CGAL_IO_READER_HELPERS_H

Polyhedron_IO/test/Polyhedron_IO/stl2off.cpp

@@ -1,53 +1,65 @@
-#include <cassert>
-#include <CGAL/IO/STL_reader.h>
-#include <CGAL/array.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
+#include <CGAL/array.h>
+#include <CGAL/IO/STL_reader.h>
+
+#include <cassert>
 #include <fstream>
 #include <iostream>
 #include <vector>
 
+template <typename Point_type, typename Polygon_type>
 void read(const char* fname, std::size_t v, std::size_t f)
 {
   std::cout << "Reading "<< fname << std::endl;
   std::ifstream input(fname, std::ios::in | std::ios::binary);
 
-  std::vector< CGAL::cpp11::array<double,3> > points;
-  std::vector< CGAL::cpp11::array<int,3> > faces;
+  std::cout << "Types: " << std::endl;
+  std::cout << typeid(Point_type).name() << std::endl;
+  std::cout << typeid(Polygon_type).name() << std::endl;
 
-  CGAL::read_STL( input,
-                  points,
-                  faces,
-                  true);
+  std::cout << "Expecting " << v << " vertices and " << f << " triangles" << std::endl;
 
-  assert(points.size() == v);
-  assert(faces.size() == f);
+  std::vector<Point_type> points;
+  std::vector<Polygon_type> faces;
+  assert(CGAL::read_STL(input, points, faces, true));
 
   std::cout << "OFF version of file " << fname << std::endl;
 
-  std::cout.precision(17);
   std::cout << "OFF\n" << points.size() << " " << faces.size()  << " 0" << std::endl;
-  for(std::size_t i=0; i < points.size(); i++){
+  for(std::size_t i=0; i < points.size(); i++)
     std::cout << points[i][0] << " " << points[i][1] << " " << points[i][2]<< std::endl;
-  }
 
-  for(std::size_t i=0; i < faces.size(); i++){
+  for(std::size_t i=0; i < faces.size(); i++)
     std::cout << "3 " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
-  }
+
+  assert(points.size() == v);
+  assert(faces.size() == f);
 }
 
-
-int main()
+int main(int, char**)
 {
-  read("data/cube.stl", 8, 12);
-  read("data/triangle.stl", 3, 1);
+  std::cout.precision(17);
 
-  read("data/ascii-tetrahedron.stl", 4, 4);
-  read("data/binary-tetrahedron-nice-header.stl", 4, 4);
-  read("data/binary-tetrahedron-non-standard-header-1.stl", 4, 4);
-  read("data/binary-tetrahedron-non-standard-header-2.stl", 4, 4);
-  read("data/binary-tetrahedron-non-standard-header-3.stl", 4, 4);
-  read("data/binary-tetrahedron-non-standard-header-4.stl", 4, 4);
-  read("data/binary-tetrahedron-non-standard-header-5.stl", 4, 4);
+  // bunch of types to test
+  typedef CGAL::cpp11::array<double, 3>                               Point_type_1;
+  typedef CGAL::Exact_predicates_exact_constructions_kernel::Point_3  Point_type_2;
+  typedef std::basic_string<double>                                   Point_type_3;
 
-  return 0;
+  typedef CGAL::cpp11::array<int, 3>                                  Polygon_type_1;
+  typedef std::vector<int>                                            Polygon_type_2;
+  typedef std::basic_string<int>                                      Polygon_type_3;
+
+  read<Point_type_1, Polygon_type_1>("data/cube.stl", 8, 12);
+  read<Point_type_1, Polygon_type_2>("data/triangle.stl", 3, 1);
+
+  read<Point_type_1, Polygon_type_3>("data/ascii-tetrahedron.stl", 4, 4);
+  read<Point_type_2, Polygon_type_1>("data/binary-tetrahedron-nice-header.stl", 4, 4);
+  read<Point_type_2, Polygon_type_2>("data/binary-tetrahedron-non-standard-header-1.stl", 4, 4);
+  read<Point_type_2, Polygon_type_3>("data/binary-tetrahedron-non-standard-header-2.stl", 4, 4);
+  read<Point_type_3, Polygon_type_1>("data/binary-tetrahedron-non-standard-header-3.stl", 4, 4);
+  read<Point_type_3, Polygon_type_2>("data/binary-tetrahedron-non-standard-header-4.stl", 4, 4);
+  read<Point_type_3, Polygon_type_3>("data/binary-tetrahedron-non-standard-header-5.stl", 4, 4);
+
+  return EXIT_SUCCESS;
 }

STL_Extension/include/CGAL/Has_member.h

@@ -0,0 +1,45 @@
+// Copyright (c) 2017 Inria (France)
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation; either version 3 of the License,
+// or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0+
+//
+// Author(s)     : Clement Jamin
+
+#ifndef CGAL_HAS_MEMBER_H
+#define CGAL_HAS_MEMBER_H
+
+// Macro used to check if a type T has a member named `X`
+// It generates a class has_X<T> where has_X<T>::value is a boolean
+// See example in Concurrent_compact_container.h
+#define CGAL_GENERATE_MEMBER_DETECTOR(X)                                      \
+template<typename T> class has_##X {                                          \
+    struct Fallback { int X; };                                               \
+    struct Derived : T, Fallback { };                                         \
+                                                                              \
+    template<typename U, U> struct Check;                                     \
+                                                                              \
+    typedef char ArrayOfOne[1];                                               \
+    typedef char ArrayOfTwo[2];                                               \
+                                                                              \
+    template<typename U> static ArrayOfOne & func(                            \
+                                            Check<int Fallback::*, &U::X> *); \
+    template<typename U> static ArrayOfTwo & func(...);                       \
+  public:                                                                     \
+    typedef has_##X type;                                                     \
+    enum { value = sizeof(func<Derived>(0)) == 2 };                           \
+} // semicolon is after the macro call
+
+#endif // CGAL_HAS_MEMBER_H

Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h

@@ -73,7 +73,7 @@ namespace CGAL {
     public:
     typedef boost::uint32_t size_type;
         /// Constructor. %Default construction creates an invalid index.
-        /// We write -1, which is <a href="http://en.cppreference.com/w/cpp/concept/numeric_limits">
+        /// We write -1, which is <a href="https://en.cppreference.com/w/cpp/types/numeric_limits">
         /// <tt>std::numeric_limits<size_type>::max()</tt></a>
         /// as `size_type` is an unsigned type. 
         explicit SM_Index(size_type _idx=(std::numeric_limits<size_type>::max)()) : idx_(_idx) {}

Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbifold_Tutte_parameterizer_3.h

@@ -210,7 +210,7 @@ Error_code read_cones(const TriangleMesh& tm, const char* filename, ConeOutputIt
 ///                  having to pass the multiple template parameters of the class `CGAL::Seam_mesh`.
 /// \tparam ConeInputBidirectionalIterator must be a model of `BidirectionalIterator`
 ///                  with value type `boost::graph_traits<SeamMesh::TriangleMesh>::%vertex_descriptor`.
-/// \tparam ConeMap must be a model of <a href="http://en.cppreference.com/w/cpp/concept/AssociativeContainer"><tt>AssociativeContainer</tt></a>
+/// \tparam ConeMap must be a model of `AssociativeContainer`
 ///                 with `boost::graph_traits<SeamMesh>::%vertex_descriptor` as key type and
 ///                 \link PkgSurfaceMeshParameterizationEnums Cone_type \endlink as value type.
 ///
@@ -896,7 +896,7 @@ public:
   /// The mapping is piecewise linear (linear in each triangle).
   /// The result is the (u,v) pair image of each vertex of the 3D surface.
   ///
-  /// \tparam ConeMap must be a model of <a href="http://en.cppreference.com/w/cpp/concept/AssociativeContainer"><tt>AssociativeContainer</tt></a>
+  /// \tparam ConeMap must be a model of `AssociativeContainer`
   ///                 with key type `boost::graph_traits<Seam_mesh>::%vertex_descriptor` and
   ///                 \link PkgSurfaceMeshParameterizationEnums Cone_type \endlink as value type.
   /// \tparam VertexUVmap must be a model of `ReadWritePropertyMap` with

Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/orbifold_shortest_path.h

@@ -175,8 +175,8 @@ void compute_shortest_paths_between_two_cones(const TriangleMesh& mesh,
 /// \tparam TriangleMesh A triangle mesh, model of `FaceListGraph` and `HalfedgeListGraph`.
 /// \tparam InputConesForwardIterator A model of `ForwardIterator` with value type
 ///                                   `boost::graph_traits<TriangleMesh>::%vertex_descriptor`.
-/// \tparam SeamContainer A model of <a href="http://en.cppreference.com/w/cpp/concept/SequenceContainer"><tt>SequenceContainer</tt></a>
-///                       with value type `boost::graph_traits<TriangleMesh>::%edge_descriptor`.
+/// \tparam SeamContainer A model of `SequenceContainer` with value type
+///                       `boost::graph_traits<TriangleMesh>::%edge_descriptor`.
 ///
 /// \param mesh the triangular mesh on which paths are computed
 /// \param first, beyond a range of cones