added hidden named parameters for isovalue nudging and interior vertex placement

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h

@@ -106,7 +106,8 @@ std::size_t get_cell_corners(const Domain& domain,
                              const typename Domain::cell_descriptor& cell,
                              const typename Domain::Geom_traits::FT isovalue,
                              Corners& corners,
-                             Values& values)
+                             Values& values,
+                             bool isovalue_nudging)
 {
   using vertex_descriptor = typename Domain::vertex_descriptor;
 
@@ -118,10 +119,10 @@ std::size_t get_cell_corners(const Domain& domain,
   for(const vertex_descriptor& v : vertices)
   {
     auto val = domain.value(v);
-    // Avoiding singular cases.
 
-    if (abs(val - isovalue) < 0.00000001)
+    // Avoiding singular cases.
-      val = isovalue - 0.00000001;
+    if (isovalue_nudging && abs(val - isovalue) < 0.00000001)
+        val = isovalue - 0.00000001;
 
     values[v_id] = val;
     if(values[v_id] >= isovalue)
@@ -130,7 +131,7 @@ std::size_t get_cell_corners(const Domain& domain,
     ++v_id;
   }
 
-  if(index.all() || index.none()) // nothing's happening in this cell
+  if(index.all() || index.none()) // nothing is happening in this cell
     return static_cast<std::size_t>(index.to_ullong());
 
   v_id = 0;
@@ -270,15 +271,18 @@ public:
 private:
   const Domain& m_domain;
   const FT m_isovalue;
+  const bool m_isovalue_nudging;
 
   Triangles m_triangles;
 
 public:
   // creates a Marching Cubes functor for a domain and isovalue
   Marching_cubes_3(const Domain& domain,
-                   const FT isovalue)
+                   const FT isovalue,
+                   const bool isovalue_nudging = true)
     : m_domain(domain),
-      m_isovalue(isovalue)
+      m_isovalue(isovalue),
+      m_isovalue_nudging(isovalue_nudging)
   { }
 
   // returns the created triangle list
@@ -300,7 +304,7 @@ public:
 
     std::array<FT, vpc> values;
     std::array<Point_3, vpc> corners;
-    const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
+    const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values, m_isovalue_nudging);
 
     // skip empty / full cells
     constexpr std::size_t ones = (1 << vpc) - 1;

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h

@@ -115,6 +115,8 @@ private:
 private:
   const Domain& m_domain;
   FT m_isovalue;
+  bool m_isovalue_nudging;
+  bool m_constrain_to_cell;
 
 #ifdef CGAL_LINKED_WITH_TBB
   std::atomic<Point_index> m_point_counter;
@@ -135,16 +137,20 @@ private:
 
 public:
   TMC_functor(const Domain& domain,
-              const FT isovalue)
+              const FT isovalue,
+              const bool isovalue_nudging = true,
+              const bool constrain_to_cell = true)
     : m_domain(domain),
       m_isovalue(isovalue),
-      m_point_counter(0)
+      m_point_counter(0),
+      m_isovalue_nudging(isovalue_nudging),
+      m_constrain_to_cell(constrain_to_cell)
   { }
 
   void operator()(const cell_descriptor& cell) {
     std::array<FT, 8> values;
     std::array<Point_3, 8> corners;
-    const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
+    const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values, m_isovalue_nudging);
 
     // skip empty / full cells
     constexpr std::size_t ones = (1 << 8) - 1;
@@ -799,7 +805,7 @@ private:
       c_ |= (size << 4 * cnt);
     };
 
-    // set corresponging edge
+    // set corresponding edge
     auto set_c = [](const int cnt, const int pos, const int val, unsigned long long& c_)
     {
       const unsigned int mask[4] = {0x0, 0xF, 0xFF, 0xFFF};
@@ -871,6 +877,7 @@ private:
       if (!std::isfinite(CGAL::to_double(ui[0])) || !std::isfinite(CGAL::to_double(ui[1])))
         continue;
 
+      // In practice, there does not seem to be a difference in choosing the first working over the best one.
       break;
     }
 
@@ -998,9 +1005,15 @@ private:
       Point_index tg_idx[6];
       for(int i=0; i<6; ++i)
       {
-        const FT u = hvt[i][0] < 0.05 ? 0.05 : (hvt[i][0] > 0.95 ? 0.95 : hvt[i][0]);
+        FT u = hvt[i][0];
-        const FT v = hvt[i][1] < 0.05 ? 0.05 : (hvt[i][1] > 0.95 ? 0.95 : hvt[i][1]);
+        FT v = hvt[i][1];
-        const FT w = hvt[i][2] < 0.05 ? 0.05 : (hvt[i][2] > 0.95 ? 0.95 : hvt[i][2]);
+        FT w = hvt[i][2];
+
+        if (m_constrain_to_cell) {
+          u = u < 0.05 ? 0.05 : (u > 0.95 ? 0.95 : u);
+          v = v < 0.05 ? 0.05 : (v > 0.95 ? 0.95 : v);
+          w = w < 0.05 ? 0.05 : (w > 0.95 ? 0.95 : w);
+        }
 
         const FT px = (FT(1) - w) * ((FT(1) - v) * (x_coord(corners[0]) + u * (x_coord(corners[1]) - x_coord(corners[0]))) +
                                                v * (x_coord(corners[2]) + u * (x_coord(corners[3]) - x_coord(corners[2])))) +
@@ -1280,9 +1293,11 @@ private:
           break;
         } // switch(c_faces)
 
-        ucoord = ucoord < 0.05 ? 0.05 : (ucoord > 0.95 ? 0.95 : ucoord);
+        if (m_constrain_to_cell) {
-        vcoord = vcoord < 0.05 ? 0.05 : (vcoord > 0.95 ? 0.95 : vcoord);
+          ucoord = ucoord < 0.05 ? 0.05 : (ucoord > 0.95 ? 0.95 : ucoord);
-        wcoord = wcoord < 0.05 ? 0.05 : (wcoord > 0.95 ? 0.95 : wcoord);
+          vcoord = vcoord < 0.05 ? 0.05 : (vcoord > 0.95 ? 0.95 : vcoord);
+          wcoord = wcoord < 0.05 ? 0.05 : (wcoord > 0.95 ? 0.95 : wcoord);
+        }
 
         // create inner vertex
         const FT px = (FT(1) - wcoord) * ((FT(1) - vcoord) * (x_coord(corners[0]) + ucoord * (x_coord(corners[1]) - x_coord(corners[0]))) +

Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h

@@ -51,6 +51,23 @@ namespace Isosurfacing {
  *     \cgalParamType{Boolean}
  *     \cgalParamDefault{`true`}
  *   \cgalParamNEnd
+ *
+ * \cond SKIP_IN_MANUAL
+ *   \cgalParamNBegin{isovalue_nudging}
+ *     \cgalParamDescription{snapping of function value at corner points to an epsilon below the isovalue if the function value is within an epsilon range around the isovalue}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`true`}
+ *     \cgalParamExtra{Snapping the function value at corner points prevents singular cases. The geometry is changed only in case of the function value being within a small epsilon range. The potential change is small. }
+ *   \cgalParamNEnd
+ *
+ *   \cgalParamNBegin{constrain_to_cell}
+ *     \cgalParamDescription{whether to exclude the cell boundary from the interior vertex placement. }
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`true`}
+ *     \cgalParamExtra{Prevents degenerate or duplicate triangles.}
+ *   \cgalParamNEnd
+ * \endcond
+ *
  * \cgalNamedParamsEnd
  *
  * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
@@ -70,17 +87,19 @@ void marching_cubes(const Domain& domain,
   using parameters::get_parameter;
 
   const bool use_tmc = choose_parameter(get_parameter(np, internal_np::use_topologically_correct_marching_cubes), true);
+  const bool isovalue_nudging = choose_parameter(get_parameter(np, internal_np::isovalue_nudging), true);
+  const bool constrain_to_cell = choose_parameter(get_parameter(np, internal_np::constrain_to_cell), true);
 
   if(use_tmc)
   {
-    internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue);
+    internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue, isovalue_nudging, constrain_to_cell);
     domain.template for_each_cell<ConcurrencyTag>(functor);
     functor.to_triangle_soup(points, triangles);
   }
   else
   {
     // run marching cubes
-    internal::Marching_cubes_3<Domain> functor(domain, isovalue);
+    internal::Marching_cubes_3<Domain> functor(domain, isovalue, isovalue_nudging);
     domain.template for_each_cell<ConcurrencyTag>(functor);
 
     // copy the result to points and triangles

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

@@ -243,6 +243,7 @@ CGAL_add_named_parameter(pca_plane_t, pca_plane, pca_plane)
 // List of named parameters used in Isosurfacing_3
 CGAL_add_named_parameter(use_topologically_correct_marching_cubes_t, use_topologically_correct_marching_cubes, use_topologically_correct_marching_cubes)
 CGAL_add_named_parameter(constrain_to_cell_t, constrain_to_cell, constrain_to_cell)
+CGAL_add_named_parameter(isovalue_nudging_t, isovalue_nudging, isovalue_nudging)
 
 // tetrahedral remeshing parameters
 CGAL_add_named_parameter(remesh_boundaries_t, remesh_boundaries, remesh_boundaries)