Tetrahedral remeshing - allow remeshing from a triangulation with badly oriented cells

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

@@ -477,7 +477,8 @@ template<typename C3t3, typename CellRange>
 bool collapse_keeps_orientations(const CellRange& star,
                                  const typename C3t3::Vertex_handle v_moved,
                                  const typename C3t3::Vertex_handle v_other,
-                                 const typename C3t3::Triangulation::Geom_traits::Point_3& new_pos)
+                                 const typename C3t3::Triangulation::Geom_traits::Point_3& new_pos,
+                                 const C3t3& c3t3)
 {
   typedef typename C3t3::Triangulation::Geom_traits::Point_3 Point;
 
@@ -490,109 +491,18 @@ bool collapse_keeps_orientations(const CellRange& star,
                                  point(ch->vertex(1)->point()),
                                  point(ch->vertex(2)->point()),
                                  point(ch->vertex(3)->point()) };
+
+    if(c3t3.triangulation().may_have_badly_oriented_cells() &&
+       CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
+      return true; //in this case we don't care about orientation change, locally
+
     pts[ch->index(v_moved)] = new_pos;
     if (CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
       return false;
   }
   return true;
 }
 
-template<typename C3t3>
-bool is_valid_collapse(const typename C3t3::Edge& edge,
-                       const Collapse_type& collapse_type,
-                       const typename C3t3::Triangulation::Point& new_pos,
-                       const C3t3& c3t3)
-{
-  typedef typename C3t3::Vertex_handle        Vertex_handle;
-  typedef typename C3t3::Cell_handle          Cell_handle;
-  typedef typename C3t3::Triangulation::Point Point;
-
-  const Vertex_handle v0 = edge.first->vertex(edge.second);
-  const Vertex_handle v1 = edge.first->vertex(edge.third);
-
-#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN
-  const bool in_cx = c3t3.is_in_complex(edge);
-  if (in_cx)
-  {
-    if (collapse_type == TO_MIDPOINT)
-      nb_test_midpoint++;
-    else if (collapse_type == TO_V1)
-      nb_test_v1++;
-    else
-      nb_test_v0++;
-  }
-#endif
-
-  if (collapse_type == TO_V1 || collapse_type == TO_MIDPOINT)
-  {
-    std::vector<Cell_handle> cells_to_check;
-    c3t3.triangulation().finite_incident_cells(v0,
-        std::back_inserter(cells_to_check));
-
-    for (const Cell_handle& ch : cells_to_check)
-    {
-      if (!ch->has_vertex(v1))
-      {
-        //check orientation
-        std::array<Point, 4> pts = { ch->vertex(0)->point(),
-                                       ch->vertex(1)->point(),
-                                       ch->vertex(2)->point(),
-                                       ch->vertex(3)->point()};
-        pts[ch->index(v0)] = new_pos;
-        if (CGAL::orientation(point(pts[0]), point(pts[1]), point(pts[2]), point(pts[3]))
-            != CGAL::POSITIVE)
-        {
-#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN
-          if (in_cx)
-          {
-            if (collapse_type == TO_MIDPOINT)
-              nb_orientation_midpoint++;
-            else
-              nb_orientation_v1++;
-          }
-#endif
-          return false;
-        }
-      }
-    }
-  }
-  if (collapse_type == TO_V0 || collapse_type == TO_MIDPOINT)
-  {
-    std::vector<Cell_handle> cells_to_check;
-    c3t3.triangulation().finite_incident_cells(v1,
-        std::back_inserter(cells_to_check));
-
-    for (const Cell_handle& ch : cells_to_check)
-    {
-      if (!ch->has_vertex(v0))
-      {
-        //check orientation
-        std::array<Point, 4> pts = { ch->vertex(0)->point(),
-                                       ch->vertex(1)->point(),
-                                       ch->vertex(2)->point(),
-                                       ch->vertex(3)->point() };
-        pts[ch->index(v1)] = new_pos;
-        if (CGAL::orientation(point(pts[0]), point(pts[1]), point(pts[2]), point(pts[3]))
-            != CGAL::POSITIVE)
-        {
-#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN
-          if (in_cx)
-          {
-            if (collapse_type == TO_MIDPOINT)
-              nb_orientation_midpoint++;
-            else
-              nb_orientation_v0++;
-          }
-#endif
-          return false;
-        }
-      }
-    }
-  }
-
-  return is_valid_collapse(edge, c3t3);
-}
-
 template<typename Facet, typename Vh>
 bool facet_has_edge(const Facet& f, const Vh v0, const Vh v1)
 {
@@ -1306,10 +1216,10 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
   const auto orientations_ok = [&](const Collapse_type ct, const Point& pos)
   {
     if ((ct == TO_V1 || ct == TO_MIDPOINT)
-        && !collapse_keeps_orientations<C3t3>(star_of_v0(), v0, v1, point(pos)))
+        && !collapse_keeps_orientations(star_of_v0(), v0, v1, point(pos), c3t3))
       return false;
     if ((ct == TO_V0 || ct == TO_MIDPOINT)
-        && !collapse_keeps_orientations<C3t3>(star_of_v1(), v1, v0, point(pos)))
+        && !collapse_keeps_orientations(star_of_v1(), v1, v0, point(pos), c3t3))
       return false;
     return true;
   };

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h

@@ -1988,7 +1988,8 @@ void flip_edges(C3T3& c3t3,
                                                    vertices_subdomain_indices);
 
 #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG
-   if(!debug::are_cell_orientations_valid(c3t3.triangulation()))
+    if(!c3t3.triangulation().may_have_badly_oriented_cells()
+     && !debug::are_cell_orientations_valid(c3t3.triangulation()))
       std::cerr << "ERROR in ORIENTATION" << std::endl;
 #endif
 

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h

@@ -496,16 +496,18 @@ class Tetrahedral_remeshing_smoother
     const typename Tr::Point backup = v->point(); //backup v's position
     const typename Tr::Geom_traits::Point_3 pv = point(backup);
 
-    bool valid_orientation = false;
-    bool angles_improved = true;
     double frac = 1.0;
     typename Tr::Geom_traits::Vector_3 move(pv, final_pos);
 
     const Dihedral_angle_cosine curr_max_cos = m_flip_smooth_steps
       ? max_cosine(tr, inc_cells)
       : Dihedral_angle_cosine(CGAL::ZERO, 0., 1.);//Dummy unused value
 
+    const bool check_orientation = need_to_check_orientation_after_change(tr, inc_cells);
+
     bool valid_try = true;
+    bool valid_orientation = true;
+    bool angles_improved = true;
     do
     {
       v->set_point(typename Tr::Point(pv + frac * move));
@@ -516,10 +518,7 @@ class Tetrahedral_remeshing_smoother
 
       for (const typename Tr::Cell_handle& ci : inc_cells)
       {
-        if (CGAL::POSITIVE != CGAL::orientation(point(ci->vertex(0)->point()),
-                                                point(ci->vertex(1)->point()),
-                                                point(ci->vertex(2)->point()),
-                                                point(ci->vertex(3)->point())))
+        if(check_orientation && !is_well_oriented(tr, ci))
         {
           frac = 0.5 * frac;
           valid_try = false;

Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h

@@ -40,10 +40,10 @@ namespace internal
 {
 
 template<typename C3t3>
-bool positive_orientation_after_edge_split(const typename C3t3::Edge& e,
+bool valid_orientation_after_edge_split(const typename C3t3::Edge& e,
                                            const typename C3t3::Cell_handle circ,
                                            const typename C3t3::Triangulation::Geom_traits::Point_3& steiner,
-                                           const C3t3&)
+                                           const C3t3& c3t3)
 {
   using Point = typename C3t3::Triangulation::Geom_traits::Point_3;
 
@@ -54,6 +54,13 @@ bool positive_orientation_after_edge_split(const typename C3t3::Edge& e,
                               point(circ->vertex(1)->point()),
                               point(circ->vertex(2)->point()),
                               point(circ->vertex(3)->point())};
+
+  // if the orientation is already bad before split, do not bother checking
+  // orientation after split, and return true straight away
+  if(c3t3.triangulation().may_have_badly_oriented_cells()
+    && CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE)
+    return true;
+
   // 1st half-cell
   const int i1 = circ->index(v1);
   const Point p1 = pts[i1];
@@ -103,7 +110,7 @@ construct_steiner_point(const typename C3t3::Edge& e,
     do
     {
       Cell_handle c = circ;
-      if(!positive_orientation_after_edge_split(e, c, steiner, c3t3))
+      if(!valid_orientation_after_edge_split(e, c, steiner, c3t3))
       {
         steiner_successful = false;
         break;
@@ -184,7 +191,7 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e,
     }
 
     const Cell_handle c = circ;
-    if(!positive_orientation_after_edge_split(e, c, m, c3t3))
+    if(!valid_orientation_after_edge_split(e, c, m, c3t3))
     {
       const std::optional<Point> steiner = construct_steiner_point(e, c3t3);
       if (steiner != std::nullopt)

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

@@ -287,8 +287,7 @@ Dihedral_angle_cosine cos_dihedral_angle(const typename Gt::Point_3& i,
                                          const typename Gt::Point_3& l,
                                          const Gt& gt)
 {
-  CGAL_expensive_assertion(CGAL::orientation(i, j, k, l) != CGAL::NEGATIVE);
-
+  //valid however orientation(i,j,k,l) is positive or negative
   typename Gt::Construct_cross_product_vector_3 cross_product =
     gt.construct_cross_product_vector_3_object();
   typename Gt::Compute_scalar_product_3 scalar_product =
@@ -365,7 +364,8 @@ Dihedral_angle_cosine max_cos_dihedral_angle(const Point& p,
   const Vector_3 ps = vector(p, s);
   const Vector_3 pr = vector(p, r);
 
-  //compute normals pointing outside tetrahedron
+  // compute normals pointing outside tetrahedron if orientation(p,q,r,s) is POSITIVE
+  // and inside tetrahedron otherwise
   const Vector_3 n_pqr = cross(qp, qr);
   if (CGAL::NULL_VECTOR == n_pqr)
     return Dihedral_angle_cosine(CGAL::POSITIVE, 1., 1.);
@@ -698,6 +698,28 @@ bool is_well_oriented(const Tr& tr, const typename Tr::Cell_handle ch)
                         ch->vertex(3));
 }
 
+template<typename Tr, typename CellRange>
+bool need_to_check_orientation_after_change(const Tr& tr, const CellRange& cells)
+{
+  for (const auto& cell : cells)
+  {
+    if (!is_well_oriented(tr, cell))
+    {
+      // if a cell is badly oriented, but the triangulation is known
+      // to have badly oriented cells, then no need to re-check orientation
+      // after change. It will remain bad
+      if(tr.may_have_badly_oriented_cells())
+        return false;
+      else
+        CGAL_assertion(false);//this function is called before the change,
+                              //so this point should not be reached
+    }
+  }
+
+  // all cells are positively oriented
+  return true;
+}
+
 template<typename C3T3, typename CellSelector>
 bool is_boundary(const C3T3& c3t3,
                  const typename C3T3::Facet& f,
@@ -2231,6 +2253,11 @@ void count_far_points(const C3t3& c3t3)
 template<typename Tr>
 bool are_cell_orientations_valid(const Tr& tr)
 {
+  // skip the test if we know
+  // that the input triangulation already has inverted cells
+  if(tr.may_have_badly_oriented_cells())
+    return true;
+
   typedef typename Tr::Geom_traits::Point_3 Point_3;
   typedef typename Tr::Facet                Facet;
 

Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp

@@ -7,26 +7,96 @@
 #include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation_3.h>
 #include <CGAL/tetrahedral_remeshing.h>
 
+#include <CGAL/SMDS_3/tet_soup_to_c3t3.h>
+#include <CGAL/Random.h>
+
 #include <iostream>
 #include <fstream>
 
-#include <CGAL/SMDS_3/tet_soup_to_c3t3.h>
-
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
 typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3<K> Remeshing_triangulation;
+typedef Remeshing_triangulation Tr;
 
-int main(int argc, char* argv[])
+
+int test(const std::string& filename, const bool perturb_vertices_positions = false)
 {
-  const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.1;
+  const double target_edge_length = 0.05;
 
-  Remeshing_triangulation tr;
+  std::ifstream in(filename);
 
-  std::ifstream in("data/sphere.mesh");
-  if (CGAL::SMDS_3::build_triangulation_from_file(in, tr))
+  Remeshing_triangulation tr;
+  if(CGAL::SMDS_3::build_triangulation_from_file(in, tr,
+    false /*verbose*/,
+    false /*replace_domain_0*/,
+    false /*allow_non_manifold*/,
+    true /*allow_negative_orientation*/))
     std::cout << "build triangulation ok" << std::endl;
+  else
+    return EXIT_FAILURE;
+
+  if (perturb_vertices_positions)
+  {
+    std::cout << "CGAL Random seed = " << CGAL::get_default_random().get_seed() << std::endl;
+    using P = Remeshing_triangulation::Point;
+    const int percentage = 10;// 1/percentage = fraction of perturbed vertices
+    const double max_move = 0.1 * target_edge_length;
+
+    CGAL::Random rng;
+    int vid = 0;
+    for (auto v : tr.finite_vertex_handles())
+    {
+      if(++vid % percentage != 0)
+        continue;
+
+      const P& pv = v->point();
+      v->set_point(P{pv.x() + rng.get_double(-1., 1.) * max_move,
+                     pv.y() + rng.get_double(-1., 1.) * max_move,
+                     pv.z() + rng.get_double(-1., 1.) * max_move});
+    }
+  }
+
+  auto count_negative_cells = [&](const Tr& tr) -> unsigned int
+    {
+      return std::count_if(tr.finite_cells_begin(),
+                           tr.finite_cells_end(),
+                           [](const auto& c)
+                           {
+                             return CGAL::POSITIVE != CGAL::orientation(c.vertex(0)->point(),
+                                                                        c.vertex(1)->point(),
+                                                                        c.vertex(2)->point(),
+                                                                        c.vertex(3)->point());
+                           });
+      };
+
+  unsigned int neg = count_negative_cells(tr);
+  if(neg > 0)
+    tr.may_have_badly_oriented_cells(true);
+  std::cout << neg << " finite cells are badly oriented before remeshing" << std::endl;
+
+  CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length,
+    CGAL::parameters::number_of_iterations(5)
+    .remesh_boundaries(false));
+
+  unsigned int neg2 = count_negative_cells(tr);
+  if(neg2 > 0)
+    std::cout << neg2 << " finite cells are badly oriented after remeshing" << std::endl;
+
+  in.close();
+  return EXIT_SUCCESS;
+}
+
+int main(int argc, char* argv[])
+{
+  std::cout << "Run remeshing from elephant.mesh" << std::endl;
+  const std::string file1 = CGAL::data_file_path("meshes/elephant.mesh");
+  if(test(file1) != EXIT_SUCCESS)
+    return EXIT_FAILURE;
 
-  CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length);
+  std::cout << "\n\nRun remeshing from perturbed elephant.mesh" << std::endl;
+  const std::string file2 = CGAL::data_file_path("meshes/elephant.mesh");
+  if(test(file2, true) != EXIT_SUCCESS)
+    return EXIT_FAILURE;
 
   return EXIT_SUCCESS;
 }