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 cabd33db179..45cfdaa47fe 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
@@ -46,6 +46,7 @@
 #include <CGAL/Isosurfacing_3/internal/tables.h>
 
 #include <CGAL/assertions.h>
+#include <CGAL/Container_helper.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
 #include <tbb/enumerable_thread_specific.h>
@@ -55,6 +56,7 @@
 
 #include <array>
 #include <bitset>
+#include <unordered_map>
 
 namespace CGAL {
 namespace Isosurfacing {
@@ -223,30 +225,40 @@ void triangles_to_polygon_soup(const TriangleRange& triangles,
                                PointRange& points,
                                PolygonRange& polygons)
 {
+  using Point = typename PointRange::value_type;
+  using PointIndexMap = std::unordered_map<Point, std::size_t>;
+
+  PointIndexMap point_index_map;
+  std::size_t current_index = 0;
+
 #ifdef CGAL_LINKED_WITH_TBB
-    for(const auto& triangle_list : triangles)
-    {
+  for(const auto& triangle_list : triangles)
+  {
 #else
-      const auto& triangle_list = triangles;
+    const auto& triangle_list = triangles;
 #endif
 
-      for(const auto& triangle : triangle_list)
+    for(const auto& triangle : triangle_list)
+    {
+      auto& polygon = polygons.emplace_back();
+      CGAL::internal::resize(polygon, 3);
+
+      for (int i=2; i>=0; --i)
       {
-        const std::size_t id = points.size();
+        const Point& p = triangle[i];
 
-        points.push_back(triangle[0]);
-        points.push_back(triangle[1]);
-        points.push_back(triangle[2]);
-
-        // simply use increasing indices
-        polygons.push_back({id + 2, id + 1, id + 0});
-
-        // just a safeguard against arrays of the wrong size
-        CGAL_assertion(polygons.back().size() == 3);
+        auto [it, inserted] = point_index_map.emplace(p, current_index);
+        if(inserted)
+        {
+          points.push_back(p);
+          ++current_index;
+        }
+        polygon[i] = it->second;
       }
+    }
 
 #ifdef CGAL_LINKED_WITH_TBB
-    }
+  }
 #endif
 }
 
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 00c69c6ef8d..f1a6ae1018a 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
@@ -80,9 +80,9 @@ private:
   using Edge_index = std::array<std::size_t, 4>;
 
 #ifdef CGAL_LINKED_WITH_TBB
-  tbb::enumerable_thread_specific<std::vector<std::array<Point_3, 3>>> m_triangles;
+  using Triangles = tbb::enumerable_thread_specific<std::vector<std::array<Point_3, 3> > >;
 #else
-  std::vector<std::array<Point_3, 3>> m_triangles;
+  using Triangles = std::vector<std::array<Point_3, 3> >;
 #endif
 
 private:
@@ -91,6 +91,8 @@ private:
   bool m_isovalue_nudging;
   bool m_constrain_to_cell;
 
+  Triangles m_triangles;
+
 public:
   TMC_functor(const Domain& domain,
               const FT isovalue,
@@ -102,6 +104,11 @@ public:
       m_constrain_to_cell(constrain_to_cell)
   { }
 
+  // returns the created triangle list
+  Triangles& triangles()
+  {
+    return m_triangles;
+  }
   void operator()(const cell_descriptor& cell) {
     std::array<FT, 8> values;
     std::array<Point_3, 8> corners;
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h
index 60cf46f19fc..4ca6d4f49ff 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h
@@ -94,15 +94,13 @@ void marching_cubes(const Domain& domain,
   {
     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);
+    internal::triangles_to_polygon_soup(functor.triangles(), points, triangles);
   }
   else
   {
     // run marching cubes
     internal::Marching_cubes_3<Domain> functor(domain, isovalue, isovalue_nudging);
     domain.template for_each_cell<ConcurrencyTag>(functor);
-
-    // copy the result to points and triangles
     internal::triangles_to_polygon_soup(functor.triangles(), points, triangles);
   }
 }