Merge pull request #6223 from janetournois/Tet_remeshing-flips_on_sur…

…face-jtournois

Tetrahedral remeshing - add edge flips on surfaces

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

@@ -3,6 +3,7 @@
 //#define CGAL_TETRAHEDRAL_REMESHING_DEBUG
 //#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS
 //#define CGAL_TETRAHEDRAL_REMESHING_PROFILE
+//#define CGAL_FLIP_ON_SURFACE_DISABLE_44_FLIP
 
 #include <QtCore/qglobal.h>
 
@@ -120,7 +121,7 @@ public Q_SLOTS:
         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);
+          Vertex_pair evv = CGAL::Tetrahedral_remeshing::make_vertex_pair(e);
           constraints.insert(evv);
         }
       }

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp

@@ -12,27 +12,27 @@
 #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;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
+using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
 
 #ifdef CGAL_CONCURRENT_MESH_3
-typedef CGAL::Parallel_tag Concurrency_tag;
+using Concurrency_tag = CGAL::Parallel_tag;
 #else
-typedef CGAL::Sequential_tag Concurrency_tag;
+using Concurrency_tag = CGAL::Sequential_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;
+using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
+using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
 
 // Criteria
-typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
 
 // Triangulation for Remeshing
-typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
-  typename Tr::Triangulation_data_structure> Triangulation_3;
+using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
+                             Tr::Triangulation_data_structure>;
 using Vertex_handle = Triangulation_3::Vertex_handle;
 
 using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
@@ -68,8 +68,7 @@ int main(int argc, char* argv[])
   // Mesh generation
   C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
 
-  CGAL::dump_c3t3(c3t3, "out");
-
+  // Property map of constraints
   Constraints_set constraints;
   Constraints_pmap constraints_pmap(constraints);
 

Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp

@@ -4,42 +4,52 @@
 #include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 #include <CGAL/Mesh_criteria_3.h>
 
+#include <CGAL/Surface_mesh.h>
 #include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
-#include <CGAL/make_mesh_3.h>
 
+#include <CGAL/make_mesh_3.h>
 #include <CGAL/tetrahedral_remeshing.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;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Polyhedron = CGAL::Surface_mesh<K::Point_3>;
+using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron>;
 
 #ifdef CGAL_CONCURRENT_MESH_3
-typedef CGAL::Parallel_tag Concurrency_tag;
+using Concurrency_tag = CGAL::Parallel_tag;
 #else
-typedef CGAL::Sequential_tag Concurrency_tag;
+using Concurrency_tag = CGAL::Sequential_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;
+// Triangulation
+using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
+using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
+  Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
 
 // Criteria
-typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
 
 // Triangulation for Remeshing
-typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
-  typename Tr::Triangulation_data_structure> Triangulation_3;
+using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
+                                              Tr::Triangulation_data_structure>;
+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");
+  const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/anchor.off");
+  const int nb_iter = (argc > 2) ? atoi(argv[2]) : 5;
+
   std::ifstream input(fname);
   Polyhedron polyhedron;
+
+  std::string filename(fname);
   input >> polyhedron;
   if (input.fail()) {
     std::cerr << "Error: Cannot read file " << fname << std::endl;
@@ -58,24 +68,29 @@ int main(int argc, char* argv[])
   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);
+  const double size = 0.072;
+  Mesh_criteria criteria(edge_size = size,
+    facet_angle = 25,
+    facet_size = size,
+    cell_radius_edge_ratio = 2,
+    cell_size = size);
 
   // Mesh generation
-  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
-
-  Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3));
-  //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));
+  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());
+
+  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));
+
+  // Remeshing
+  CGAL::tetrahedral_isotropic_remeshing(tr, size,
+    CGAL::parameters::number_of_iterations(nb_iter));
+
+  std::ofstream out("out_remeshed.mesh");
+  CGAL::IO::write_MEDIT(out, tr);
+  out.close();
 
   return EXIT_SUCCESS;
 }

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h

@@ -44,6 +44,7 @@ enum Edge_type     { FEATURE, BOUNDARY, INSIDE, MIXTE,
 enum Collapse_type { TO_MIDPOINT, TO_V0, TO_V1, IMPOSSIBLE };
 enum Result_type   { VALID,
                      V_PROBLEM, C_PROBLEM, E_PROBLEM,
+                     ANGLE_PROBLEM,
                      TOPOLOGICAL_PROBLEM, ORIENTATION_PROBLEM, SHARED_NEIGHBOR_PROBLEM };
 
 template<typename C3t3>
@@ -214,6 +215,15 @@ class CollapseTriangulation
 
       } while (++circ != done);
 
+#ifdef PROTECT_ANGLES_FROM_COLLAPSE
+      Dihedral_angle_cosine curr_max_cos
+        = max_cos_dihedral_angle(triangulation, cells_to_remove[0]);
+      for (std::size_t i = 1; i < cells_to_remove.size(); ++i)
+      {
+        curr_max_cos = (std::max)(curr_max_cos,
+          max_cos_dihedral_angle(triangulation, cells_to_remove[i]));
+      }
+#endif
 
       vh0->set_point(Point_3(v0_new_pos.x(), v0_new_pos.y(), v0_new_pos.z()));
       vh1->set_point(Point_3(v0_new_pos.x(), v0_new_pos.y(), v0_new_pos.z()));
@@ -269,6 +279,10 @@ class CollapseTriangulation
           return ORIENTATION_PROBLEM;
         if (!triangulation.tds().is_valid(cit, true))
           return C_PROBLEM;
+#ifdef PROTECT_ANGLES_FROM_COLLAPSE
+        if (curr_max_cos < max_cos_dihedral_angle(triangulation, cit))
+          return ANGLE_PROBLEM;
+#endif
       }
 
       for (Vertex_handle vit : triangulation.finite_vertex_handles())
@@ -791,6 +805,9 @@ collapse(const typename C3t3::Cell_handle ch,
   }
   while (++circ != done);
 
+  if(c3t3.is_in_complex(ch->vertex(from), ch->vertex(to)))
+    c3t3.remove_from_complex(ch->vertex(from), ch->vertex(to));
+
   bool valid = true;
   std::vector<Cell_handle> cells_to_remove;
   std::unordered_set<Cell_handle> invalid_cells;
@@ -1109,11 +1126,9 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
   if (are_edge_lengths_valid(edge, c3t3, new_pos, sqhigh, cell_selector/*, adaptive = false*/)
     && collapse_preserves_surface_star(edge, c3t3, new_pos, cell_selector))
   {
-    CGAL_assertion_code(typename Tr::Cell_handle dc);
-    CGAL_assertion_code(int di);
-    CGAL_assertion_code(int dj);
-    CGAL_assertion(c3t3.triangulation().is_edge(edge.first->vertex(edge.second),
-                                                edge.first->vertex(edge.third), dc, di, dj));
+    CGAL_assertion(c3t3.triangulation().tds().is_edge(
+                       edge.first->vertex(edge.second),
+                       edge.first->vertex(edge.third)));
 
     Vertex_handle v0_init = edge.first->vertex(edge.second);
     Vertex_handle v1_init = edge.first->vertex(edge.third);
@@ -1223,7 +1238,7 @@ void collapse_short_edges(C3T3& c3t3,
 
     FT sqlen = sql(tr.segment(e));
     if (sqlen < sq_low)
-      short_edges.insert(short_edge(make_vertex_pair<T3>(e), sqlen));
+      short_edges.insert(short_edge(make_vertex_pair(e), sqlen));
   }
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
@@ -1285,7 +1300,7 @@ void collapse_short_edges(C3T3& c3t3,
 
           const FT sqlen = sql(tr.segment(eshort));
           if (sqlen < sq_low)
-            short_edges.insert(short_edge(make_vertex_pair<T3>(eshort), sqlen));
+            short_edges.insert(short_edge(make_vertex_pair(eshort), sqlen));
         }
 
         //debug::dump_c3t3(c3t3, "dump_after_collapse");