Merge duplicate points in (T)MC post processing

Pointless to leave that step to the user.

Maybe the soup builder could do that on-the-fly, but then you have
to use concurrent data structures, not sure that it would be faster?

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,6 +225,12 @@ 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)
   {
@@ -232,17 +240,21 @@ void triangles_to_polygon_soup(const TriangleRange& triangles,
 
     for(const auto& triangle : triangle_list)
     {
-        const std::size_t id = points.size();
+      auto& polygon = polygons.emplace_back();
+      CGAL::internal::resize(polygon, 3);
 
-        points.push_back(triangle[0]);
+      for (int i=2; i>=0; --i)
-        points.push_back(triangle[1]);
+      {
-        points.push_back(triangle[2]);
+        const Point& p = triangle[i];
 
-        // simply use increasing indices
+        auto [it, inserted] = point_index_map.emplace(p, current_index);
-        polygons.push_back({id + 2, id + 1, id + 0});
+        if(inserted)
-
+        {
-        // just a safeguard against arrays of the wrong size
+          points.push_back(p);
-        CGAL_assertion(polygons.back().size() == 3);
+          ++current_index;
+        }
+        polygon[i] = it->second;
+      }
     }
 
 #ifdef CGAL_LINKED_WITH_TBB

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;

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);
   }
 }