fix property map for region primitives

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh_planar_patches.h

@@ -165,13 +165,14 @@ struct Face_index_tracker<TriangleMeshOut, VertexCornerMapOut, internal_np::Para
   Fmap f2f_map() { return Fmap(); }
 };
 
-template <class Vector_3, class NormalRange>
+template <class Vector_3, class PatchNormalMap>
 void init_face_normals(std::vector<Vector_3>& face_normals,
                        std::size_t nb_patches,
-                       const NormalRange& normal_range)
+                       PatchNormalMap patch_normal_map)
 {
-  face_normals.reserve(nb_patches);
+  face_normals.resize(nb_patches);
-  face_normals.assign(normal_range.begin(), normal_range.end());
+  for (std::size_t i=0; i<nb_patches; ++i)
+    face_normals[i] = get(patch_normal_map, i);
 }
 
 template <class Vector_3>
@@ -741,7 +742,18 @@ bool decimate_impl(const TriangleMesh& tm,
         {
           //TODO: shall we try to plug pseudo-cdt?
 #ifdef CGAL_DEBUG_DECIMATION
-          std::cout << "  DEBUG: Failed to remesh a patch" << std::endl;
+          static int fail_case_id=0;
+          std::cout << "  DEBUG: Failed to remesh a patch, case #" << fail_case_id  << std::endl;
+          std::ofstream debug("failed_remesh_"+std::to_string(fail_case_id)+".polylines.txt");
+          debug << std::setprecision(17);
+          for (auto c : csts)
+            debug << "2 " << corners[c.first] << " " << corners[c.second] << "\n";
+          debug.close();
+          std::cout << "  normal used is " << face_normals[cc_id] << "\n";
+          debug.open("normal"+std::to_string(fail_case_id)+".polylines.txt");
+          debug << "2 " << corners[csts[0].first] << " " << corners[csts[0].first]+face_normals[cc_id] << "\n";
+          debug.close();
+          ++fail_case_id;
 #endif
           all_patches_successfully_remeshed = false;
           // make all vertices of the patch a corner
@@ -1501,6 +1513,17 @@ void remesh_planar_patches(const TriangleMeshIn& tm_in,
  *  \param np_in an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below:
  *
  *  \cgalNamedParamsBegin
+ *    \cgalParamNBegin{patch_normal_map}
+ *      \cgalParamDescription{a property map providing for each patch the normal of the supporting plane of the patch (used to triangulate it)}
+ *      \cgalParamType{a class model of `ReadPropertyMap` with the value type of `FacePatchMap` as key and
+ *                     `GeomTraits::Vector_3` as value type, `GeomTraits` being the type of the parameter `geom_traits`}
+ *      \cgalParamDefault{If not provided, patch normals will be estimated using corners of the patches}
+ *    \cgalParamNEnd
+ *    \cgalParamNBegin{do_not_triangulate_faces}
+ *      \cgalParamDescription{if `true`, faces of `out` will not be triangulated, but the one with more than one connected component of the boundary.}
+ *      \cgalParamType{`bool`}
+ *      \cgalParamDefault{false}
+ *    \cgalParamNEnd
  *    \cgalParamNBegin{vertex_point_map}
  *      \cgalParamDescription{a property map associating points to the vertices of `tm_in`}
  *      \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<TriangleMeshIn>::%vertex_descriptor`
@@ -1509,11 +1532,6 @@ void remesh_planar_patches(const TriangleMeshIn& tm_in,
  *      \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
  *                      must be available in `TriangleMeshIn`.}
  *    \cgalParamNEnd
- *    \cgalParamNBegin{do_not_triangulate_faces}
- *      \cgalParamDescription{if `true`, faces of `out` will not be triangulated, but the one with more than one connected component of the boundary.}
- *      \cgalParamType{`bool`}
- *      \cgalParamDefault{false}
- *   \cgalParamNEnd
  *    \cgalParamNBegin{geom_traits}
  *      \cgalParamDescription{an instance of a geometric traits class}
  *      \cgalParamType{a class model of `Kernel`}
@@ -1540,12 +1558,6 @@ void remesh_planar_patches(const TriangleMeshIn& tm_in,
  *                     as key type and `std::size_t` as value type}
  *      \cgalParamDefault{None}
  *    \cgalParamNEnd
- *    \cgalParamNBegin{normals_of_patches}
- *      \cgalParamDescription{the normals of the supporting planes of the patches (in the same order of `face_patch_map`) used to triangulate patches.}
- *      \cgalParamType{a model of `InputIterator` with `GeomTraits::Vector_3` as value type,
- *                     `GeomTraits` being the type of the parameter `geom_traits`}
- *      \cgalParamDefault{normals will be estimated from three non-collinear corners on the boundary of each patch}
- *    \cgalParamNEnd
  *    \cgalParamNBegin{vertex_corner_map}
  *      \cgalParamDescription{a property map filled by this function and that will contain for each vertex its corner
  *                            an id in the range `[0, number of corners - 1]`}
@@ -1599,7 +1611,7 @@ bool remesh_almost_planar_patches(const TriangleMeshIn& tm_in,
   bool do_not_triangulate_faces = choose_parameter(get_parameter(np_in, internal_np::do_not_triangulate_faces), false);
 
   std::vector< typename Traits::Vector_3 > face_normals;
-  Planar_segmentation::init_face_normals(face_normals, nb_patches, get_parameter(np_out, internal_np::normals_of_patches));
+  Planar_segmentation::init_face_normals(face_normals, nb_patches, get_parameter(np_in, internal_np::patch_normal_map));
   return Planar_segmentation::decimate_impl<Traits>(tm_in, tm_out,
                                                     std::make_pair(nb_corners, nb_patches),
                                                     vertex_corner_map, ecm, face_patch_map, vpm_in, vpm_out,

Polyhedron/demo/Polyhedron/Plugins/PMP/Remesh_planar_patches_plugin.cpp

@@ -13,6 +13,7 @@
 #include <CGAL/property_map.h>
 
 #include <boost/graph/graph_traits.hpp>
+#include <boost/property_map/vector_property_map.hpp>
 
 #include <QElapsedTimer>
 #include <QAction>
@@ -198,12 +199,12 @@ public Q_SLOTS:
         Patch_id_pmap in_fpmap = get(CGAL::face_patch_id_t<int>(), pmesh);
         std::vector<std::size_t> corner_id_map(num_vertices(pmesh), -1);
         std::vector<bool> ecm(num_edges(pmesh), false);
-        std::vector<CGAL::Epick::Vector_3> normals;
+        boost::vector_property_map<CGAL::Epick::Vector_3> normal_map;
         std::size_t nb_regions =
           PMP::region_growing_of_planes_on_faces(pmesh,
                                                  in_fpmap,
-                                                 normals,
                                                  CGAL::parameters::cosine_of_maximum_angle(cos_threshold).
+                                                                   region_primitive_map(normal_map).
                                                                    maximum_distance(ui.dist_dspinbox->value()));
         std::size_t nb_corners =
           PMP::detect_corners_of_regions(pmesh,
@@ -227,7 +228,7 @@ public Q_SLOTS:
                                             CGAL::make_random_access_property_map(corner_id_map),
                                             CGAL::make_random_access_property_map(ecm),
                                             CGAL::parameters::do_not_triangulate_faces(do_not_triangulate_faces).
-                                                              normals_of_patches(CGAL::make_random_access_property_map(normals)),
+                                                              patch_normal_map(normal_map),
                                             CGAL::parameters::face_patch_map(out_fpmap));
 
 
@@ -253,7 +254,7 @@ public Q_SLOTS:
                                             CGAL::make_random_access_property_map(corner_id_map),
                                             CGAL::make_random_access_property_map(ecm),
                                             CGAL::parameters::do_not_triangulate_faces(do_not_triangulate_faces).
-                                                              normals_of_patches(CGAL::make_random_access_property_map(normals)),
+                                                              patch_normal_map(normal_map),
                                             CGAL::parameters::visitor([](Mesh& pmesh){pmesh.clear_without_removing_property_maps ();}));
           pmesh.remove_property_map<Mesh::Face_index, int>(in_fpmap);
           poly_item->invalidateOpenGLBuffers();

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -148,7 +148,8 @@ CGAL_add_named_parameter(mesh_facet_distance_t, mesh_facet_distance, mesh_facet_
 CGAL_add_named_parameter(mesh_facet_topology_t, mesh_facet_topology, mesh_facet_topology)
 CGAL_add_named_parameter(polyline_constraints_t, polyline_constraints, polyline_constraints)
 CGAL_add_named_parameter(vertex_corner_map_t, vertex_corner_map, vertex_corner_map)
-CGAL_add_named_parameter(normals_of_patches_t, normals_of_patches, normals_of_patches)
+CGAL_add_named_parameter(patch_normal_map_t, patch_normal_map, patch_normal_map)
+CGAL_add_named_parameter(region_primitive_map_t, region_primitive_map, region_primitive_map)
 
 // 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)

Shape_detection/include/CGAL/Polygon_mesh_processing/region_growing.h

@@ -25,6 +25,32 @@
 namespace CGAL {
 namespace Polygon_mesh_processing {
 
+namespace internal
+{
+  template <class GT, class Pair>
+  void fill_region_primitive_map(const std::vector<Pair>&, internal_np::Param_not_found) {}
+
+  template <class GT, class Pair, class RegionMap>
+  void fill_plane_or_vector_map(const std::vector<Pair>& normals, RegionMap region_map, typename GT::Vector_3)
+  {
+    for (std::size_t i = 0 ; i<normals.size(); ++i)
+      put(region_map, i, normals[i].first.orthogonal_vector());
+  }
+
+  template <class GT, class Pair, class RegionMap>
+  void fill_plane_or_vector_map(const std::vector<Pair>& normals, RegionMap region_map, typename GT::Plane_3)
+  {
+    for (std::size_t i = 0 ; i<normals.size(); ++i)
+      put(region_map, i, normals[i].first);
+  }
+
+  template <class GT, class Pair, class RegionMap>
+  void fill_region_primitive_map(const std::vector<Pair>& normals, RegionMap region_map)
+  {
+    fill_plane_or_vector_map<GT>(normals, region_map, typename boost::property_traits<RegionMap>::value_type());
+  }
+}
+
 /*!
   \ingroup PkgPolygonMeshProcessingRef
   \brief applies a region growing algorithm to fit planes on faces of a mesh.
@@ -81,6 +107,14 @@ namespace Polygon_mesh_processing {
       \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
       \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
     \cgalParamNEnd
+    \cgalParamNBegin{region_primitive_map}
+      \cgalParamDescription{a property map filled by this function and that will contain for each region
+                            the plane (or only its orthognonal vector) estimated that approximates it.}
+      \cgalParamType{a class model of `WritablePropertyMap` with the value type of `RegionMap` as key and
+                     `GeomTraits::Plane_3` or `GeomTraits::Vector_3` as value type,
+                     `GeomTraits` being the type of the parameter `geom_traits`}
+      \cgalParamDefault{None}
+    \cgalParamNEnd
   \cgalNamedParamsEnd
  */
 template<typename PolygonMesh,
@@ -89,7 +123,6 @@ template<typename PolygonMesh,
 std::size_t
 region_growing_of_planes_on_faces(const PolygonMesh& mesh,
                                   RegionMap region_map,
-                                  std::vector<typename GetGeomTraits<PolygonMesh, NamedParameters>::type::Vector_3> & normals,
                                   const NamedParameters& np = parameters::default_values())
 {
   namespace RG_PM = CGAL::Shape_detection::Polygon_mesh;
@@ -118,11 +151,7 @@ region_growing_of_planes_on_faces(const PolygonMesh& mesh,
   std::vector<std::pair<typename Traits::Plane_3, std::vector<typename PolygonMesh::Face_index> > > tmp;
   region_growing.detect(std::back_inserter(tmp));
 
-  normals.resize(region_growing.number_of_regions_detected());
+  internal::fill_region_primitive_map<Traits>(tmp, parameters::get_parameter(np, internal_np::region_primitive_map));
-  for (std::size_t i =0 ; i<region_growing.number_of_regions_detected(); ++i)
-  {
-    normals[i]=tmp[i].first.orthogonal_vector();
-  }
 
   return region_growing.number_of_regions_detected();
 }