Merge branch 'gsoc2022-isosurface' of https://github.com/JulyCode/cgal into gsoc2022-isosurface

# Conflicts:
#	Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp

@@ -5,6 +5,8 @@
 #include <CGAL/boost/graph/IO/OFF.h>
 #include <vector>
 
+#include <tbb/tick_count.h>
+
 using Kernel = CGAL::Simple_cartesian<double>;
 using FT = typename Kernel::FT;
 using Point = typename Kernel::Point_3;
@@ -32,10 +34,16 @@ int main(int, char**)
   Point_range points;
   Triangle_range triangles;
 
+  const tbb::tick_count start = tbb::tick_count::now();
+
   // execute marching cubes with a given isovalue 
   const FT isovalue = 0.8;
   CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles);
 
+  const tbb::tick_count end = tbb::tick_count::now();
+
+  std::cout << (end - start).seconds() << std::endl;
+
   // save ouput indexed mesh to a file, in the OFF format
   CGAL::IO::write_OFF("output.off", points, triangles);
 

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp

@@ -17,7 +17,7 @@ using Polygon_range = std::vector<std::vector<std::size_t> >;
 
 int main(int, char**)
 {
-  const std::string fname = CGAL::data_file_path("images/skull_2.9.inr");
+  const std::string fname = "../examples/Isosurfacing_3/FullHead.inr";//CGAL::data_file_path("images/skull_2.9.inr");
 
   // load volumetric image from a file
   CGAL::Image_3 image;
@@ -30,6 +30,11 @@ int main(int, char**)
   // convert image to a Cartesian grid
   Grid grid{image};
 
+  for (std::size_t i = 0; i < grid.xdim(); i++)
+    for (std::size_t j = 0; j < grid.ydim(); j++)
+      for (std::size_t k = 0; k < grid.zdim(); k++)
+        grid.value(i, j, k) = 2 * 1120 - grid.value(i, j, k);
+
   // create a domain from the grid
   auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid);
 
@@ -38,7 +43,7 @@ int main(int, char**)
   Polygon_range polygons;
 
   // execute marching cubes
-  CGAL::Isosurfacing::marching_cubes(domain, 2.9 /*isovalue*/, points, polygons);
+  CGAL::Isosurfacing::marching_cubes(domain, 1120 /*isovalue*/, points, polygons);
 
   // save output indexed mesh to a file, in the OFF format
   CGAL::IO::write_OFF("result.off", points, polygons);

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h

@@ -20,6 +20,7 @@
 
 #ifdef CGAL_LINKED_WITH_TBB
 #include <tbb/parallel_for.h>
+#include <tbb/blocked_range3d.h>
 #endif // CGAL_LINKED_WITH_TBB
 
 #include <array>
@@ -213,14 +214,22 @@ public:
     const std::size_t sk = size_k;
 
     // for now only parallelize outer loop
-    auto iterator = [&f, sj, sk](const tbb::blocked_range<std::size_t>& r) {
+    auto iterator = [&f, sj, sk](const tbb::blocked_range3d<std::size_t>& r) {
-      for(std::size_t i = r.begin(); i != r.end(); ++i)
+      const std::size_t i_begin = r.pages().begin();
-        for(std::size_t j=0; j<sj - 1; ++j)
+      const std::size_t i_end = r.pages().end();
-          for(std::size_t k=0; k<sk - 1; ++k)
+      const std::size_t j_begin = r.rows().begin();
+      const std::size_t j_end = r.rows().end();
+      const std::size_t k_begin = r.cols().begin();
+      const std::size_t k_end = r.cols().end();
+
+      for(std::size_t i = i_begin; i != i_end; ++i)
+        for(std::size_t j = j_begin; j != j_end; ++j)
+          for(std::size_t k = k_begin; k != k_end; ++k)
             f({i, j, k});
     };
 
-    tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_i - 1), iterator);
+    tbb::blocked_range3d<std::size_t> range(0, size_i - 1, 0, size_j - 1, 0, size_k - 1);
+    tbb::parallel_for(range, iterator);
   }
 #endif // CGAL_LINKED_WITH_TBB
 

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

@@ -48,7 +48,7 @@
 #include <CGAL/assertions.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
-#include <tbb/concurrent_vector.h>
+#include <tbb/enumerable_thread_specific.h>
 #else
 #include <vector>
 #endif
@@ -187,7 +187,12 @@ void mc_construct_triangles(const int i_case,
     const int eg2 = Cube_table::triangle_cases[t_index + 2];
 
     // insert new triangle in list
-    triangles.push_back({vertices[eg0], vertices[eg1], vertices[eg2]});
+    #ifdef CGAL_LINKED_WITH_TBB
+      auto& tris = triangles.local();
+    #else
+      auto& tris = triangles;
+    #endif
+    tris.push_back({vertices[eg0], vertices[eg1], vertices[eg2]});
   }
 }
 
@@ -198,17 +203,28 @@ void triangles_to_polygon_soup(const TriangleRange& triangles,
                                PointRange& points,
                                PolygonRange& polygons)
 {
-  for(auto& triangle : triangles)
+  #ifdef CGAL_LINKED_WITH_TBB
-  {
+    for(const auto& triangle_list : triangles)
-    const std::size_t id = points.size();
+    {
+  #else
+      const auto& triangle_list = triangles;
+  #endif
 
-    points.push_back(triangle[0]);
+      for(const auto& triangle : triangle_list)
-    points.push_back(triangle[1]);
+      {
-    points.push_back(triangle[2]);
+        const std::size_t id = points.size();
 
-    // simply use increasing indices
+        points.push_back(triangle[0]);
-    polygons.push_back({id + 2, id + 1, id + 0});
+        points.push_back(triangle[1]);
-  }
+        points.push_back(triangle[2]);
+
+        // simply use increasing indices
+        polygons.push_back({id + 2, id + 1, id + 0});
+      }
+
+  #ifdef CGAL_LINKED_WITH_TBB
+    }
+  #endif
 }
 
 // Marching Cubes implemented as a functor that runs on every cell of the grid
@@ -225,7 +241,7 @@ public:
   using Cell_descriptor = typename Domain::Cell_descriptor;
 
 #ifdef CGAL_LINKED_WITH_TBB
-  using Triangles = tbb::concurrent_vector<std::array<Point_3, 3> >;
+  using Triangles = tbb::enumerable_thread_specific<std::vector<std::array<Point_3, 3>>>;
 #else
   using Triangles = std::vector<std::array<Point_3, 3> >;
 #endif
@@ -245,7 +261,7 @@ public:
   { }
 
   // gets the created triangle list
-  const Triangles& triangles() const
+  Triangles& triangles()
   {
     return m_triangles;
   }

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

@@ -46,11 +46,16 @@
 #include <CGAL/Isosurfacing_3/internal/marching_cubes_functors.h>
 #include <CGAL/Isosurfacing_3/internal/tables.h>
 
+#include <tbb/concurrent_vector.h>
+#include <tbb/concurrent_hash_map.h>
+
 #include <array>
 #include <atomic>
 #include <cmath>
 #include <map>
 #include <mutex>
+#include <atomic>
+#include <optional>
 
 namespace CGAL {
 namespace Isosurfacing {
@@ -73,26 +78,45 @@ private:
   using Edge_descriptor = typename Domain::Edge_descriptor;
   using Cell_descriptor = typename Domain::Cell_descriptor;
 
-  using uint = unsigned int;
+  using Point_index = std::size_t;
+  using Edge_index = std::array<std::size_t, 4>;
+
+  struct Hash_compare
+  {
+    static size_t hash(const Edge_index& key)
+    {
+      std::size_t res = 17;
+      res = res * 31 + std::hash<std::size_t>()(key[0]);
+      res = res * 31 + std::hash<std::size_t>()(key[1]);
+      res = res * 31 + std::hash<std::size_t>()(key[2]);
+      res = res * 31 + std::hash<std::size_t>()(key[3]);
+      return res;
+    }
+
+    static bool equal(const Edge_index& key1, const Edge_index& key2)
+    {
+      return key1[0] == key2[0] && key1[1] == key2[1] && key1[2] == key2[2] && key1[3] == key2[3];
+    }
+  };
+
+  using Edge_point_map = tbb::concurrent_hash_map<Edge_index, Point_index, Hash_compare>;
 
 private:
   const Domain& m_domain;
   FT m_isovalue;
 
-  Point_range& m_points;
+  std::atomic<Point_index> m_point_counter;
-  Polygon_range& m_polygons;
+  tbb::concurrent_vector<Point_3> m_points;
 
-  std::mutex mutex;
+  Edge_point_map m_edges;
+
+  tbb::concurrent_vector<std::array<Point_index, 3>> m_triangles;
 
 public:
   TMC_functor(const Domain& domain,
-              const FT isovalue,
+              const FT isovalue)
-              Point_range& points,
-              Polygon_range& polygons)
     : m_domain(domain),
-      m_isovalue(isovalue),
+      m_isovalue(isovalue)
-      m_points(points),
-      m_polygons(polygons)
   { }
 
   void operator()(const Cell_descriptor& cell)
@@ -101,18 +125,19 @@ public:
     std::array<Point_3, 8> corners;
     const int i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values);
 
-    const int all_bits_set = (1 << (8 + 1)) - 1;  // last 8 bits are 1
-    if(Cube_table::intersected_edges[i_case] == 0 ||
-       Cube_table::intersected_edges[i_case] == all_bits_set)
-    {
-      return;
-    }
-
     // this is the only difference to mc
-    int tcm = int(Cube_table::t_ambig[i_case]);
+    const int tcm = Cube_table::t_ambig[i_case];
     if(tcm == 105)
     {
-      p_slice(cell, m_isovalue, values, corners, i_case);
+      if (p_slice(cell, m_isovalue, values, corners, i_case))
+        return;
+      else
+        std::cerr << "WARNING: the result might not be topologically correct" << std::endl;
+    }
+
+    constexpr int all_bits_set = (1 << (8 + 1)) - 1;  // last 8 bits are 1
+    if(i_case == 0 || i_case == all_bits_set)
+    {
       return;
     }
 
@@ -122,7 +147,6 @@ public:
     // @todo improve triangle generation
 
     // construct triangles
-    std::lock_guard<std::mutex> lock(mutex);
     for(int t=0; t<16; t += 3)
     {
       const int t_index = i_case * 16 + t;
@@ -135,38 +159,90 @@ public:
       const int eg1 = Cube_table::triangle_cases[t_index + 1];
       const int eg2 = Cube_table::triangle_cases[t_index + 2];
 
-      const std::size_t p0_idx = m_points.size();
-
-      m_points.push_back(vertices[eg0]);
-      m_points.push_back(vertices[eg1]);
-      m_points.push_back(vertices[eg2]);
 
       // insert new triangle into list
-      m_polygons.emplace_back();
+      const Point_index p0 = add_point(vertices[eg0], compute_edge_index(cell, eg0));
-      auto& triangle = m_polygons.back();
+      const Point_index p1 = add_point(vertices[eg1], compute_edge_index(cell, eg1));
+      const Point_index p2 = add_point(vertices[eg2], compute_edge_index(cell, eg2));
 
-      triangle.push_back(p0_idx + 2);
+      add_triangle(p2, p1, p0);
-      triangle.push_back(p0_idx + 1);
+    }
-      triangle.push_back(p0_idx + 0);
+  }
+
+  // gets the created triangle list
+  template<typename PointRange, typename TriangleRange>
+  void to_triangle_soup(PointRange& points, TriangleRange& triangles) const
+  {
+    points.insert(points.begin(), m_points.begin(), m_points.end());
+    for (const std::array<Point_index, 3>& tri : m_triangles) {
+      triangles.push_back({ tri[0], tri[1], tri[2] });
     }
   }
 
 private:
-  void add_triangle(const std::size_t p0,
+  Edge_index compute_edge_index(const Cell_descriptor& cell, int edge)
-                    const std::size_t p1,
-                    const std::size_t p2)
   {
-    std::lock_guard<std::mutex> lock(mutex);
+    // edge is in 0 - 11
 
-    m_polygons.emplace_back();
+    // there are 12 edges, assign to each vertex three edges, the global edge numbering
-    auto& triangle = m_polygons.back();
+    // consists of 3*global_vertex_id + edge_offset.
+    const unsigned long long gei_pattern_ = 670526590282893600ull;
 
-    triangle.push_back(p0);
+    // the edge global index is given by the vertex global index + the edge offset
-    triangle.push_back(p1);
+    const std::size_t shift = 5 * edge;
-    triangle.push_back(p2);
+    const std::size_t ix = cell[0] + ((gei_pattern_ >> shift) & 1);        // global_edge_id[edge][0];
+    const std::size_t iy = cell[1] + ((gei_pattern_ >> (shift + 1)) & 1);  // global_edge_id[edge][1];
+    const std::size_t iz = cell[2] + ((gei_pattern_ >> (shift + 2)) & 1);  // global_edge_id[edge][2];
+    const std::size_t off_val =      ((gei_pattern_ >> (shift + 3)) & 3);
+
+    return { ix, iy, iz, off_val };
   }
 
-  void p_slice(const Cell_descriptor& cell,
+  bool find_point(const Edge_index& e, Point_index& i)
+  {
+    Edge_point_map::const_accessor acc;
+    if (m_edges.find(acc, e))
+    {
+      i = acc->second;
+      return true;
+    }
+    return false;
+  }
+
+  Point_index add_point(const Point_3& p, const Edge_index& e)
+  {
+    Edge_point_map::accessor acc;
+    if (!m_edges.insert(acc, e))
+      return acc->second;
+    
+    const Point_index i = m_point_counter++;
+    acc->second = i;
+    acc.release();
+
+    m_points.grow_to_at_least(i + 1);
+    m_points[i] = p;
+
+    return i;
+  }
+
+  Point_index add_point_unchecked(const Point_3& p)
+  {
+    const Point_index i = m_point_counter++;
+
+    m_points.grow_to_at_least(i + 1);
+    m_points[i] = p;
+
+    return i;
+  }
+
+  void add_triangle(const Point_index p0,
+                    const Point_index p1,
+                    const Point_index p2)
+  {
+    m_triangles.push_back({p0, p1, p2});
+  }
+
+  bool p_slice(const Cell_descriptor& cell,
                const FT i0,
                const std::array<FT, 8>& values,
                const std::array<Point_3, 8>& corners,
@@ -177,9 +253,7 @@ private:
     typename Geom_traits::Compute_z_3 z_coord = m_domain.geom_traits().compute_z_3_object();
     typename Geom_traits::Construct_point_3 point = m_domain.geom_traits().construct_point_3_object();
 
-    // there are 12 edges, assign to each vertex three edges, the global edge numbering
+    using uint = unsigned int;
-    // consists of 3*global_vertex_id + edge_offset.
-    const unsigned long long gei_pattern_ = 670526590282893600ull;
 
     // code edge end vertices for each of the 12 edges
     const unsigned char l_edges_[12] = {16, 49, 50, 32, 84, 117, 118, 100, 64, 81, 115, 98};
@@ -191,25 +265,18 @@ private:
 
     // A hexahedron has twelve edges, save the intersection of the isosurface with the edge
     // save global edge and global vertex index of isosurface
-    std::vector<std::size_t> vertices(12);
+    std::vector<Point_index> vertices(12);
 
     // save local coordinate along the edge of intersection point
     std::vector<FT> ecoord{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
     // collect vertices
     unsigned short flag{1};
-    for(int eg=0; eg<12; ++eg)
+    for(int eg = 0; eg < 12; ++eg)
     {
       if(flag & Cube_table::intersected_edges[i_case])
       {
-        // the edge global index is given by the vertex global index + the edge offset
+        
-        // uint shift = 5 * eg;
-        // const int ix = i_index + (int)((gei_pattern_ >> shift) & 1);        // global_edge_id[eg][0];
-        // const int iy = j_index + (int)((gei_pattern_ >> (shift + 1)) & 1);  // global_edge_id[eg][1];
-        // const int iz = k_index + (int)((gei_pattern_ >> (shift + 2)) & 1);  // global_edge_id[eg][2];
-        // const int off_val = (int)((gei_pattern_ >> (shift + 3)) & 3);
-
-        // int g_edg = int(m_cell_shift_factor * m_ugrid.global_index(ix, iy, iz) + off_val);
 
         // generate vertex here, do not care at this point if vertex already exists
         uint v0, v1;
@@ -223,22 +290,9 @@ private:
         const FT py = (1 - l) * y_coord(corners[v0]) + l * y_coord(corners[v1]);
         const FT pz = (1 - l) * z_coord(corners[v0]) + l * z_coord(corners[v1]);
 
-        // set vertex in map
+        // add vertex and insert to map
-        // set vertex index
+        vertices[eg] = add_point(point(px, py, pz), compute_edge_index(cell, eg));
-        // auto s_index = m_vertices.find(vertices[eg].g_edg);
-        // if(s_index == m_vertices.end())
-        // {
-          const int g_idx = static_cast<int>(m_points.size());
-          vertices[eg] = g_idx;
-        //   m_vertices[vertices[eg].g_edg] = g_idx;
-          m_points.push_back(point(px, py, pz));
-        //} else {
-        //   vertices[eg] = s_index->second;
-        //}
       }
-      /*else {
-        e_set[eg] = false;
-      }*/
 
       // next edge
       flag <<= 1;
@@ -409,8 +463,8 @@ private:
             }
             else
             {
-              std::cerr << "ERROR: can't correctly triangulate cell's face\n";
+              // std::cerr << "ERROR: can't correctly triangulate cell's face\n";
-              return;
+              return false;
             }
           }
         }
@@ -481,8 +535,8 @@ private:
             }
              else
             {
-              std::cerr << "ERROR: can't correctly triangulate cell's face\n";
+              // std::cerr << "ERROR: can't correctly triangulate cell's face\n";
-              return;
+              return false;
             }
           }
         }
@@ -815,7 +869,7 @@ private:
 
       // compute vertices of inner hexagon, save new vertices in list and compute and keep
       // global vertices index to build triangle connectivity later on.
-      std::size_t tg_idx[6];
+      Point_index tg_idx[6];
       for(int i=0; i<6; ++i)
       {
         const FT u = hvt[i][0];
@@ -834,8 +888,7 @@ private:
                             w * ((1 - v) * (z_coord(corners[4]) + u * (z_coord(corners[5]) - z_coord(corners[4]))) +
                                        v * (z_coord(corners[6]) + u * (z_coord(corners[7]) - z_coord(corners[6]))));
 
-        tg_idx[i] = m_points.size();
+        tg_idx[i] = add_point_unchecked(point(px, py, pz));
-        m_points.push_back(point(px, py, pz));
       }
 
       // triangulate contours with inner hexagon
@@ -1113,10 +1166,10 @@ private:
                             wcoord * ((1 - vcoord) * (z_coord(corners[4]) + ucoord * (z_coord(corners[5]) - z_coord(corners[4]))) +
                                             vcoord * (z_coord(corners[6]) + ucoord * (z_coord(corners[7]) - z_coord(corners[6]))));
 
-        const std::size_t g_index = m_points.size();
+        bool pt_used = false;
+        Point_index g_index = 0;
 
         // loop over the contours
-        bool pt_used = false;
         for(int i=0; i<int(cnt_); ++i)
         {
           const unsigned char cnt_sz = (unsigned char)get_cnt_size(i, c_);
@@ -1126,18 +1179,21 @@ private:
           }
           else
           {
-            pt_used = true;
+            if (!pt_used)
+            {
+              pt_used = true;
+              g_index = add_point_unchecked(point(px, py, pz));
+            }
+            
             for(int t=0; t<cnt_sz; ++t)
             {
               add_triangle(vertices[get_c(i, t, c_)], vertices[get_c(i, (t + 1) % cnt_sz, c_)], g_index);
             }
           }
         }
-
-        if(pt_used)
-          m_points.emplace_back(px, py, pz);
       }
     }
+    return true;
   }
 };
 

Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h

@@ -74,8 +74,9 @@ void marching_cubes(const Domain& domain,
   {
     // run topologically correct marching cubes
     // and directly write the result to points and triangles
-    internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue, points, triangles);
+    internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue);
     domain.template iterate_cells<ConcurrencyTag>(functor);
+    functor.to_triangle_soup(points, triangles);
   }
   else
   {