Tetrahedral remeshing - deal with c3t3 complex edges (#7909)

## Summary of Changes

For a mesh generated by Mesh_3, that has complex edges that do not lie
at the intersection of surface patches, automatic detection of these
features is impossible from the input triangulation.
This PR introduces a new named parameter `edge_is_constrained_map` to
`convert_to_triangulation_3()` that sets the edges constrained status
corresponding to their `is_in_complex()` status.

This pmap can then be used straight away in
`tetrahedral_isotropic_remeshing()`.

## Release Management

* Affected package(s): Tetrahedral_remeshing
* License and copyright ownership: unchanged

Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp

@@ -12,6 +12,7 @@
 #include "C3t3_type.h"
 
 #include <CGAL/tetrahedral_remeshing.h>
+#include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h>
 
 #include <unordered_map>
 #include <memory>
@@ -80,6 +81,7 @@ public Q_SLOTS:
 
     Scene_c3t3_item* c3t3_item =
       qobject_cast<Scene_c3t3_item*>(scene->item(index));
+    const auto& c3t3 = c3t3_item->c3t3();
 
     if (c3t3_item)
     {
@@ -100,18 +102,37 @@ public Q_SLOTS:
       bool protect = ui.protect_checkbox->isChecked();
       bool smooth_edges = ui.smoothEdges_checkBox->isChecked();
 
+      // collect constraints
+      using Vertex_handle = Tr::Vertex_handle;
+      using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
+      using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
+      using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
+
       // wait cursor
       QApplication::setOverrideCursor(Qt::WaitCursor);
 
       QElapsedTimer time;
       time.start();
 
+      Constraints_set constraints;
+      for (const auto e : c3t3_item->c3t3().triangulation().finite_edges())
+      {
+        if (c3t3_item->c3t3().is_in_complex(e)
+          || CGAL::Tetrahedral_remeshing::protecting_balls_intersect(e, c3t3))
+        {
+          Vertex_pair evv = CGAL::Tetrahedral_remeshing::make_vertex_pair<Tr>(e);
+          constraints.insert(evv);
+        }
+      }
+
       CGAL::tetrahedral_isotropic_remeshing(
-          c3t3_item->c3t3(),
-          target_length,
-          CGAL::parameters::remesh_boundaries(!protect)
-          .number_of_iterations(nb_iter)
-          .smooth_constrained_edges(smooth_edges));
+        c3t3_item->c3t3(),
+        target_length,
+        CGAL::parameters::remesh_boundaries(!protect)
+        .number_of_iterations(nb_iter)
+        .smooth_constrained_edges(smooth_edges)
+        .edge_is_constrained_map(Constraints_pmap(constraints))
+      );
 
       std::cout << "Remeshing done (" << time.elapsed() << " ms)" << std::endl;
 

Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt

@@ -103,6 +103,12 @@ setting the named parameter `remesh_boundaries` to `false`.
 
 \cgalExample{Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp }
 
+It is also possible to define the polyline features as the ones
+stored as complex edges in a `Mesh_complex_3_in_triangulation_3`
+(e.g. generated by the \ref PkgMesh3 package mesh generation algorithms).
+
+\cgalExample{Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp }
+
 
 \subsection ssecEx4 Tetrahedral Remeshing After Mesh Generation
 

Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt

@@ -4,5 +4,6 @@
 \example Tetrahedral_remeshing/tetrahedral_remeshing_of_one_subdomain.cpp
 \example Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp
 \example Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp
+\example Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp
 \example Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp
 */

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt

@@ -31,9 +31,13 @@ if(TARGET CGAL::Eigen3_support)
   create_single_source_cgal_program( "mesh_and_remesh_polyhedral_domain_with_features.cpp" )
   target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen3_support)
 
+  create_single_source_cgal_program("mesh_and_remesh_c3t3.cpp")
+  target_link_libraries(mesh_and_remesh_c3t3 PUBLIC CGAL::Eigen3_support)
+
   if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
     message(STATUS "Found TBB")
     target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PRIVATE CGAL::TBB_support)
+    target_link_libraries(mesh_and_remesh_c3t3 PRIVATE CGAL::TBB_support)
   endif()
 else()
   message(STATUS "NOTICE: Some examples require Eigen 3.1 (or greater), and will not be compiled.")

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp

@@ -0,0 +1,96 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+
+#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
+#include <CGAL/make_mesh_3.h>
+
+#include <CGAL/tetrahedral_remeshing.h>
+
+#include <CGAL/IO/File_medit.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
+
+#ifdef CGAL_CONCURRENT_MESH_3
+typedef CGAL::Parallel_tag Concurrency_tag;
+#else
+typedef CGAL::Sequential_tag Concurrency_tag;
+#endif
+
+// Triangulation for Meshing
+typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;
+typedef CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+// Triangulation for Remeshing
+typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
+  typename Tr::Triangulation_data_structure> Triangulation_3;
+using Vertex_handle = Triangulation_3::Vertex_handle;
+
+using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
+using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
+using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
+
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main(int argc, char* argv[])
+{
+  const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fandisk.off");
+  std::ifstream input(fname);
+  Polyhedron polyhedron;
+  input >> polyhedron;
+  if (input.fail() || !CGAL::is_triangle_mesh(polyhedron)) {
+    std::cerr << "Error: Input invalid " << fname << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // Create domain
+  Mesh_domain domain(polyhedron);
+
+  // Get sharp features
+  domain.detect_features();
+
+  // Mesh criteria
+  Mesh_criteria criteria(edge_size = 0.025,
+    facet_angle = 25, facet_size = 0.05, facet_distance = 0.005,
+    cell_radius_edge_ratio = 3, cell_size = 0.05);
+
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
+
+  CGAL::dump_c3t3(c3t3, "out");
+
+  Constraints_set constraints;
+  Constraints_pmap constraints_pmap(constraints);
+
+  Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3),
+    CGAL::parameters::edge_is_constrained_map(constraints_pmap));
+
+  //note we use the move semantic, with std::move(c3t3),
+  //  to avoid a copy of the triangulation by the function
+  //  `CGAL::convert_to_triangulation_3()`
+  //  After the call to this function, c3t3 is an empty and valid C3t3.
+  //It is possible to use :  CGAL::convert_to_triangulation_3(c3t3),
+  //  Then the triangulation is copied and duplicated, and c3t3 remains as is.
+
+  const double target_edge_length = 0.1;//coarsen the mesh
+  CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length,
+    CGAL::parameters::number_of_iterations(3)
+   .edge_is_constrained_map(constraints_pmap));
+
+  std::ofstream out("out_remeshed.mesh");
+  CGAL::IO::write_MEDIT(out, tr);
+  out.close();
+
+  return EXIT_SUCCESS;
+}

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp

@@ -12,59 +12,16 @@
 
 #include "tetrahedral_remeshing_generate_input.h"
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
 
-typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
+using Remeshing_triangulation = CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K>;
+using Point = Remeshing_triangulation::Point;
+using Vertex_handle = Remeshing_triangulation::Vertex_handle;
 
-typedef Remeshing_triangulation::Point         Point;
-typedef Remeshing_triangulation::Vertex_handle Vertex_handle;
-typedef Remeshing_triangulation::Cell_handle   Cell_handle;
-typedef Remeshing_triangulation::Edge          Edge;
+using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
+using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
+using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
 
-class Constrained_edges_property_map
-{
-public:
-  typedef bool                               value_type;
-  typedef bool                               reference;
-  typedef std::pair<Vertex_handle, Vertex_handle> key_type;
-  typedef boost::read_write_property_map_tag category;
-
-private:
-  std::unordered_set<key_type, boost::hash<key_type>>* m_set_ptr;
-
-public:
-  Constrained_edges_property_map()
-    : m_set_ptr(nullptr)
-  {}
-  Constrained_edges_property_map(std::unordered_set<key_type, boost::hash<key_type>>* set_)
-    : m_set_ptr(set_)
-  {}
-
-public:
-  friend void put(Constrained_edges_property_map& map,
-                  const key_type& k,
-                  const bool b)
-  {
-    assert(map.m_set_ptr != nullptr);
-    assert(k.first < k.second);
-    if (b)  map.m_set_ptr->insert(k);
-    else    map.m_set_ptr->erase(k);
-  }
-
-  friend value_type get(const Constrained_edges_property_map& map,
-                        const key_type& k)
-  {
-    assert(map.m_set_ptr != nullptr);
-    assert(k.first < k.second);
-    return (map.m_set_ptr->count(k) > 0);
-  }
-};
-
-void set_subdomain(Remeshing_triangulation& tr, const int index)
-{
-  for (auto cit : tr.finite_cell_handles())
-    cit->set_subdomain_index(index);
-}
 
 int main(int argc, char* argv[])
 {
@@ -73,16 +30,14 @@ int main(int argc, char* argv[])
   const int nbv = (argc > 3) ? atoi(argv[3]) : 500;
 
   Remeshing_triangulation t3;
-  typedef std::pair<Vertex_handle, Vertex_handle> Vertex_pair;
-  std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>> constraints;
 
+  Constraints_set constraints;
   CGAL::Tetrahedral_remeshing::generate_input_cube(nbv, t3, constraints);
   make_constraints_from_cube_edges(t3, constraints);
 
   CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length,
-    CGAL::parameters::edge_is_constrained_map(
-      Constrained_edges_property_map(&constraints))
-    .number_of_iterations(nb_iter));
+    CGAL::parameters::edge_is_constrained_map(Constraints_pmap(constraints))
+                     .number_of_iterations(nb_iter));
 
   return EXIT_SUCCESS;
 }

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h

@@ -471,17 +471,8 @@ private:
     const Curve_index default_curve_id = default_curve_index();
     for (const Edge& e : tr().finite_edges())
     {
-      if (m_c3t3.is_in_complex(e))
-      {
-        CGAL_assertion(m_c3t3.in_dimension(e.first->vertex(e.second)) <= 1);
-        CGAL_assertion(m_c3t3.in_dimension(e.first->vertex(e.third)) <= 1);
-#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-        ++nbe;
-#endif
-        continue;
-      }
-
       if (get(ecmap, CGAL::Tetrahedral_remeshing::make_vertex_pair<Tr>(e))
+          || m_c3t3.is_in_complex(e)
           || nb_incident_subdomains(e, m_c3t3) > 2
           || nb_incident_surface_patches(e, m_c3t3) > 1)
       {
@@ -539,6 +530,7 @@ private:
     CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension(
       m_c3t3.triangulation(), "0-c3t3_vertices_after_init_");
     CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3);
+    CGAL::Tetrahedral_remeshing::debug::count_far_points(m_c3t3);
 #endif
   }
 

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h

@@ -568,6 +568,8 @@ void set_index(typename C3t3::Vertex_handle v, const C3t3& c3t3)
   case 0:
     v->set_index(Mesh_3::internal::get_index<typename C3t3::Corner_index>(v->index()));
     break;
+  case -1://far points from concurrent Mesh_3
+    break;
   default:
     CGAL_assertion(false);
   }
@@ -589,6 +591,27 @@ bool is_edge_in_complex(const typename C3t3::Vertex_handle& v0,
     return false;
 }
 
+template<typename C3t3>
+bool protecting_balls_intersect(const typename C3t3::Edge& e,
+                                const C3t3& c3t3)
+{
+  const auto vv = c3t3.triangulation().vertices(e);
+  if(  c3t3.in_dimension(vv[0]) > 1
+    || c3t3.in_dimension(vv[1]) > 1)
+    return false;
+
+  const auto& p0 = vv[0]->point();
+  const auto& p1 = vv[1]->point();
+  if(p0.weight() == 0 || p1.weight() == 0)
+    return false;
+
+  const auto r0 = CGAL::approximate_sqrt(p0.weight());
+  const auto r1 = CGAL::approximate_sqrt(p1.weight());
+  const auto d = CGAL::approximate_sqrt(CGAL::squared_distance(p0, p1));
+
+  return d < r0 + r1;
+}
+
 template<typename C3t3, typename OutputIterator>
 OutputIterator incident_subdomains(const typename C3t3::Vertex_handle v,
                                    const C3t3& c3t3,
@@ -1268,6 +1291,18 @@ void check_surface_patch_indices(const C3t3& c3t3)
   }
 }
 
+template<typename C3t3>
+void count_far_points(const C3t3& c3t3)
+{
+  std::size_t count = 0;
+  for (auto v : c3t3.triangulation().finite_vertex_handles())
+  {
+    if(c3t3.in_dimension(v) == -1)
+      ++count;
+  }
+  std::cout << "Nb far points : " << count << std::endl;
+}
+
 template<typename Tr>
 bool are_cell_orientations_valid(const Tr& tr)
 {
@@ -1678,13 +1713,17 @@ void dump_vertices_by_dimension(const Tr& tr, const char* prefix)
   typedef typename Tr::Vertex_handle Vertex_handle;
   std::vector< std::vector<Vertex_handle> > vertices_per_dimension(4);
 
+  std::size_t nb_far_points = 0;
   for (typename Tr::Finite_vertices_iterator
        vit = tr.finite_vertices_begin();
        vit != tr.finite_vertices_end();
        ++vit)
   {
     if (vit->in_dimension() == -1)
+    {
+      ++nb_far_points;
       continue;//far point
+    }
     CGAL_assertion(vit->in_dimension() >= 0 && vit->in_dimension() < 4);
 
     vertices_per_dimension[vit->in_dimension()].push_back(vit);
@@ -1712,6 +1751,7 @@ void dump_vertices_by_dimension(const Tr& tr, const char* prefix)
 
     ofs.close();
   }
+  std::cout << "Nb far points : " << nb_far_points << std::endl;
 }
 
 template<typename Tr>

Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h

@@ -310,21 +310,64 @@ void tetrahedral_isotropic_remeshing(
 * @tparam CurveIndex is the type of the indices for feature curves.
 *            If `c3t3` has been generated using `CGAL::make_mesh_3()`, it must match
 *            `MeshDomainWithFeatures_3::Curve_index`.
+* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
 *
 * @param c3t3 the complex containing the triangulation to be remeshed.
+* @param np optional sequence of \ref bgl_namedparameters "Named Parameters"
+*          among the ones listed below
+*
+* \cgalNamedParamsBegin
+*   \cgalParamNBegin{edge_is_constrained_map}
+*     \cgalParamDescription{a property map containing the constrained-or-not status of each edge of
+*                     `c3t3.triangulation()`.
+*                     For each edge `e` for which `c3t3.is_in_complex(e)` returns `true`,
+*                     the constrained status of `e` is set to `true`.}
+*     \cgalParamType{a class model of `ReadWritePropertyMap`
+*         with `std::pair<Triangulation_3::Vertex_handle, Triangulation_3::Vertex_handle>`
+*         as key type and `bool` as value type. It must be default constructible.}
+*     \cgalParamDefault{a default property map where no edge is constrained}
+*   \cgalParamNEnd
+*
+* \cgalNamedParamsEnd
 */
 
 template<typename Tr,
          typename CornerIndex,
-         typename CurveIndex>
+         typename CurveIndex,
+         typename NamedParameters = parameters::Default_named_parameters>
 CGAL::Triangulation_3<typename Tr::Geom_traits,
                       typename Tr::Triangulation_data_structure>
 convert_to_triangulation_3(
-  CGAL::Mesh_complex_3_in_triangulation_3<Tr, CornerIndex, CurveIndex> c3t3)
+  CGAL::Mesh_complex_3_in_triangulation_3<Tr, CornerIndex, CurveIndex> c3t3,
+  const NamedParameters& np = parameters::default_values())
 {
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
   using GT   = typename Tr::Geom_traits;
   using TDS  = typename Tr::Triangulation_data_structure;
 
+  using Vertex_handle = typename Tr::Vertex_handle;
+  using Edge_vv = std::pair<Vertex_handle, Vertex_handle>;
+  using Default_pmap = Constant_property_map<Edge_vv, bool>;
+  using ECMap = typename internal_np::Lookup_named_param_def <
+                internal_np::edge_is_constrained_t,
+                NamedParameters,
+                Default_pmap
+                >::type;
+  ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
+                                 Default_pmap(false));
+
+  if (!std::is_same_v<ECMap, Default_pmap>)
+  {
+    for (auto e : c3t3.edges_in_complex())
+    {
+      const Edge_vv evv
+        = CGAL::Tetrahedral_remeshing::make_vertex_pair<Tr>(e);//ordered pair
+      put(ecmap, evv, true);
+    }
+  }
+
   CGAL::Triangulation_3<GT, TDS> tr;
   tr.swap(c3t3.triangulation());
   return tr;