Avoid computing points when the cell is trivial

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

@@ -97,7 +97,8 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi
                z_coord(p1) * mu + z_coord(p0) * (FT(1) - mu));
 }
 
-// retrieves the corner vertices and their values of a cell and return the lookup index
+// retrieves the values of a cell and return the lookup index
+// if the cell is completely above or below the isovalue, corner points are not computed
 template <typename Domain,
           typename Corners,
           typename Values>
@@ -109,22 +110,27 @@ std::size_t get_cell_corners(const Domain& domain,
 {
   using Vertex_descriptor = typename Domain::Vertex_descriptor;
 
+  const auto& vertices = domain.cell_vertices(cell);
+
   // collect function values and build index
   std::size_t v_id = 0;
   std::bitset<Domain::VERTICES_PER_CELL> index = 0;
-  for(const Vertex_descriptor& v : domain.cell_vertices(cell))
+  for(const Vertex_descriptor& v : vertices)
   {
-    // collect scalar values and computex index
-    corners[v_id] = domain.point(v);
     values[v_id] = domain.value(v);
-
     if(values[v_id] >= isovalue)
       index.set(v_id);
 
-    // next cell vertex
     ++v_id;
   }
 
+  if(index.all() || index.none()) // nothing's happening in this cell
+    return static_cast<std::size_t>(index.to_ullong());
+
+  v_id = 0;
+  for(const Vertex_descriptor& v : vertices)
+    corners[v_id++] = domain.point(v);
+
   return static_cast<std::size_t>(index.to_ullong());
 }
 
@@ -290,7 +296,7 @@ public:
     std::array<Point_3, vpc> corners;
     const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
 
-    // skip empty cells
+    // skip empty / full cells
     constexpr std::size_t ones = (1 << vpc) - 1;
     if((i_case & ones) == ones || // all bits set
        (i_case & ones) == 0) // no bits set

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

@@ -146,11 +146,17 @@ public:
     std::array<Point_3, 8> corners;
     const std::size_t i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
 
+    // skip empty / full cells
+    constexpr std::size_t ones = (1 << 8) - 1;
+    if((i_case & ones) == ones || // all bits set
+       (i_case & ones) == 0) // no bits set
+      return;
+
     // this is the only difference to the default Marching Cubes
     const int tcm = Cube_table::t_ambig[i_case];
     if(tcm == 105)
     {
-      if(p_slice(cell, m_isovalue, values, corners, i_case))
+      if(p_slice(cell, m_isovalue, corners, values, i_case))
         return;
 #ifdef CGAL_ISOSURFACING_3_MC_FUNCTORS_DEBUG
       else
@@ -158,11 +164,6 @@ public:
 #endif
     }
 
-    constexpr std::size_t ones = (1 << 8) - 1;
-    if((i_case & ones) == ones || // all bits set
-       (i_case & ones) == 0) // no bits set
-      return;
-
     std::array<Point_3, 12> vertices;
     MC_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices);
 
@@ -294,8 +295,8 @@ private:
 
   bool p_slice(const Cell_descriptor& cell,
                const FT i0,
-               const std::array<FT, 8>& values,
                const std::array<Point_3, 8>& corners,
+               const std::array<FT, 8>& values,
                const int i_case)
   {
     typename Geom_traits::Compute_x_3 x_coord = m_domain.geom_traits().compute_x_3_object();