diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h
index b9a7d0eac3c..c9489890bc1 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h
@@ -78,20 +78,14 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi
   typename GeomTraits::Compute_z_3 z_coord = gt.compute_z_3_object();
   typename GeomTraits::Construct_point_3 point = gt.construct_point_3_object();
 
-  FT mu = FT(0);
-
   // @todo, technically we should be using the edge intersection oracle here, but there is a nuance
   // between MC and DC on the handling of edges that have val0 = val1 = isovalue: in MC we assume
   // the isosurface is in the middle, in DC we assume the isosurface is not intersecting the edge.
   // In the oracle, we follow DC right now. Could put a Boolean parameter, but it's ugly.
 
-  // don't divide by 0
-  if(abs(d1 - d0) < 0.000001) // @fixme hardcoded bound
-    mu = FT(0.5);  // if both points have the same value, assume isolevel is in the middle
-  else
-    mu = (isovalue - d0) / (d1 - d0);
-
-  CGAL_assertion(mu >= FT(0.0) || mu <= FT(1.0));
+  const FT den = d1 - d0;
+  FT mu = is_zero(den) ? FT(1) / FT(2) : (isovalue - d0) / den;
+  mu = std::clamp<FT>(mu, FT(0), FT(1));
 
   // linear interpolation
   return point((FT(1) - mu) * x_coord(p0) +  mu * x_coord(p1),
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h
index b5d31c0621b..d4e0e30d887 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h
@@ -400,7 +400,10 @@ private:
         get_edge_vertex(eg, v0, v1, l_edges_);
 
         // @todo use the domain's interpolation scheme?
-        FT l = (i0 - values[v0]) / (values[v1] - values[v0]);
+        const FT den = values[v1] - values[v0];
+        FT l = is_zero(den) ? FT(1) / FT(2) : (i0 - values[v0]) / den;
+        l = std::clamp<FT>(l, FT(0), FT(1));
+
         ecoord[eg] = l;
 
         const FT px = (FT(1) - l) * x_coord(corners[v0]) + l * x_coord(corners[v1]);