Tetrahedral remeshing - deal with c3t3 complex edges

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 @@ class Adaptive_remesher
     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 @@ class Adaptive_remesher
     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>