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

Documentation/doc/CMakeLists.txt

@@ -131,6 +131,10 @@ function(configure_doxygen_package CGAL_PACKAGE_NAME)
       endif()
     endif()
   endif()
+  if(EXISTS "${CGAL_PACKAGE_DIR}/include/CGAL/${CGAL_PACKAGE_NAME}/internal")
+    file(APPEND ${CGAL_DOC_PACKAGE_DEFAULTS}
+        "EXCLUDE += ${CGAL_PACKAGE_DIR}/include/CGAL/${CGAL_PACKAGE_NAME}/internal\n")
+  endif()
 
   # IMAGE_PATH is set by default. For Documentation, we generate the extra path using packages.txt
   set(IMAGE_PATHS "${CGAL_PACKAGE_DOC_DIR}/fig")

Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt

@@ -41,10 +41,10 @@ if(TARGET CGAL::TBB_support)
   target_link_libraries(marching_cubes_inrimage PRIVATE CGAL::TBB_support)
 
   if(TARGET CGAL::Eigen3_support)
-	target_link_libraries(dual_contouring_cartesian_grid PRIVATE CGAL::TBB_support)
+    target_link_libraries(dual_contouring_cartesian_grid PRIVATE CGAL::TBB_support)
-	target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::TBB_support)
+    target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::TBB_support)
-	target_link_libraries(dual_contouring_octree PRIVATE CGAL::TBB_support)
+    target_link_libraries(dual_contouring_octree PRIVATE CGAL::TBB_support)
-	target_link_libraries(all_cartesian_cube PRIVATE CGAL::TBB_support)
+    target_link_libraries(all_cartesian_cube PRIVATE CGAL::TBB_support)
-	target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::TBB_support)
+    target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::TBB_support)
   endif()
 endif()

Isosurfacing_3/include/CGAL/Dual_contouring_3.h

@@ -21,6 +21,31 @@
 namespace CGAL {
 namespace Isosurfacing {
 
+/**
+ * \ingroup PkgIsosurfacing3Ref
+ *
+ * \brief Creates an indexed face set that represents an isosurface using the Dual Contouring algorithm.
+ *
+ * \details
+ *
+ * \tparam ConcurrencyTag determines if the algorithm is executed sequentially or in parallel.
+ *
+ * \tparam Domain_ must be a model of `IsosurfacingDomain`.
+ *
+ * \tparam PointRange is a model of the concept `RandomAccessContainer` and `BackInsertionSequence` whose value type can
+ * be constructed from the point type of the polygon mesh.
+ * \tparam PolygonRange a model of the concept
+ * `RandomAccessContainer` and `BackInsertionSequence` whose value type is itself a model of the concepts
+ * `RandomAccessContainer` and `BackInsertionSequence` whose value type is `std::size_t`.
+ * \tparam Positioning is a functor containing the operator() that takes `domain`, `iso_value`, `cell`, and `position`
+ * as input and returns a boolean that is true if the isosurface intersects the cell.
+ *
+ * \param domain the domain providing input data and its topology
+ * \param iso_value value of the isosurface
+ * \param points points making the polygons of the indexed face set
+ * \param polygons each element in the vector describes a polygon using the indices of the points in points
+ * \param positioning the functor dealing with vertex positioning inside a voxel
+ */
 template <typename Concurrency_tag = Sequential_tag, class Domain_, class PointRange, class PolygonRange,
           class Positioning = internal::Positioning::QEM_SVD<true>>
 void dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value, PointRange& points,
@@ -28,25 +53,29 @@ void dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value
 
     // static_assert(Domain_::CELL_TYPE & ANY_CELL);
 
-    internal::Dual_contouring_position_functor<Domain_, Positioning> pos_func(domain, iso_value, positioning);
+    internal::Dual_contouring_vertex_positioning<Domain_, Positioning> pos_func(domain, iso_value, positioning);
     domain.iterate_cells(pos_func, Concurrency_tag());
 
-    internal::Dual_contouring_quads_functor<Domain_> quad_func(domain, iso_value);
+    internal::Dual_contouring_face_generation<Domain_> face_generation(domain, iso_value);
-    domain.iterate_edges(quad_func, Concurrency_tag());
+    domain.iterate_edges(face_generation, Concurrency_tag());
 
-    // write points and quads in ranges
+    // write points and faces in ranges
     points.resize(pos_func.points_counter);
     for (const auto& vtop : pos_func.map_voxel_to_point) {
         points[pos_func.map_voxel_to_point_id[vtop.first]] = vtop.second;
     }
 
-    polygons.reserve(quad_func.quads.size());
+    polygons.reserve(face_generation.faces.size());
-    for (const auto& q : quad_func.quads) {
+    for (const auto& q : face_generation.faces) {
         std::vector<std::size_t> vertex_ids;
         for (const auto& v_id : q.second) {
-            vertex_ids.push_back(pos_func.map_voxel_to_point_id[v_id]);
+            if (pos_func.map_voxel_to_point_id.count(v_id) > 0) {
+                vertex_ids.push_back(pos_func.map_voxel_to_point_id[v_id]);
+            }
+        }
+        if (vertex_ids.size() > 2) {
+            polygons.push_back(vertex_ids);
         }
-        polygons.push_back(vertex_ids);
     }
 }
 

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Dual_contouring_internal.h

@@ -13,13 +13,13 @@
 #ifndef CGAL_DUAL_CONTOURING_3_INTERNAL_DUAL_CONTOURING_3_H
 #define CGAL_DUAL_CONTOURING_3_INTERNAL_DUAL_CONTOURING_3_H
 
-#include <CGAL/license/Isosurfacing_3.h>
 #include <CGAL/Bbox_3.h>
-#include <CGAL/Isosurfacing_3/internal/Tables.h>
 #include <CGAL/Origin.h>
+#include <CGAL/centroid.h>
+#include <CGAL/license/Isosurfacing_3.h>
 
-#include <array>
 #include <Eigen/SVD>
+#include <array>
 #include <map>
 #include <mutex>
 #include <vector>
@@ -29,121 +29,112 @@ namespace Isosurfacing {
 namespace internal {
 
 namespace Positioning {
+/**
+ * \ingroup PkgIsosurfacing3Ref
+ *
+ * \brief Computes the vertex position for a point in Dual Contouring using Quadric Error Metrics and the SVD pseudo
+ * inverse.
+ *
+ * \details
+ *
+ * \tparam use_bbox clamp vertex position to the bounding box of the cell
+ *
+ */
 template <bool use_bbox = false>
 class QEM_SVD {
 public:
-    /// <summary>
+    /**
-    /// Compute vertex position for Dual Contouring
+     * \ingroup PkgIsosurfacing3Ref
-    /// </summary>
+     *
-    /// <typeparam name="Domain_"></typeparam>
+     * \brief Computes the vertex position for a point in Dual Contouring.
-    /// <param name="domain"></param>
+     *
-    /// <param name="iso_value"></param>
+     * \details
-    /// <param name="i"></param>
+     *
-    /// <param name="j"></param>
+     * \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
-    /// <param name="k"></param>
+     *
-    /// <returns> true, if there is a point in the cell</returns>
+     * \param domain the domain providing input data and its topology
+     * \param iso_value value of the isosurface
+     * \param cell the cell within the domain for which the vertex position ins computed
+     * \param point the point position of the vertex that belongs to that cell
+     *
+     * \return true, if the voxel intersects the isosurface
+     */
     template <class Domain_>
-    bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh,
+    bool position(const Domain_& domain, const typename Domain_::FT iso_value,
-                  typename Domain_::Point& point) const {
+                  const typename Domain_::Cell_handle& cell, typename Domain_::Point& point) const {
         typedef typename Domain_::Point Point;
         typedef typename Domain_::Geom_traits::Vector_3 Vector;
         typedef typename Domain_::FT FT;
 
-        typename Domain_::Cell_vertices vertices = domain.cell_vertices(vh);
+        typename Domain_::Cell_vertices vertices = domain.cell_vertices(cell);
 
-        namespace Tables = internal::Cube_table;
+        std::vector<Point> pos(vertices.size());
-
-        std::array<typename Domain_::FT, Tables::N_VERTICES> s;
-        std::transform(vertices.begin(), vertices.end(), s.begin(), [&](const auto& v) { return domain.value(v); });
-
-        std::array<bool, Tables::N_VERTICES> b;
-        std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; });
-
-        unsigned int cubeindex = 0;
-        // set bit if corresponding corner is below iso
-        for (int i = 0; i < Tables::N_VERTICES; ++i) {
-            cubeindex |= b[i] << i;
-        }
-
-        if (cubeindex == 0 || cubeindex == 255) {
-            return false;
-        }
-
-        std::array<Vector, Tables::N_VERTICES> pos;
         std::transform(vertices.begin(), vertices.end(), pos.begin(),
-                       [&](const auto& v) { return domain.position(v) - CGAL::ORIGIN; });
+                       [&](const auto& v) { return domain.position(v); });
 
-        point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0]));  // set point to voxel center
+        point = CGAL::centroid(pos.begin(), pos.end(), CGAL::Dimension_tag<0>());  // set point to cell center
-
-        // std::array<Vector, Tables::N_VERTICES> normals;
-        // std::transform(vertices.begin(), vertices.end(), normals.begin(),
-        //                [&](const auto& v) { return domain.gradient(domain.position(v)); });
 
         // compute edge intersections
         std::vector<Point> edge_intersections;
         std::vector<Vector> edge_intersection_normals;
 
-        for (int i = 0; i < Tables::N_EDGES; ++i) {
+        for (const auto& edge : domain.cell_edges(cell)) {
-            const auto& v0 = Tables::edge_to_vertex[i][0];
+            const auto& edge_vertices = domain.edge_vertices(edge);
-            const auto& v1 = Tables::edge_to_vertex[i][1];
+            const auto& v0 = edge_vertices[0];
+            const auto& v1 = edge_vertices[1];
 
-            if (b[v0] != b[v1]) {  // e0
+            const auto& val0 = domain.value(v0);
-                const FT u = (s[v0] - iso_value) / (s[v0] - s[v1]);
+            const auto& val1 = domain.value(v1);
-                const Point p_lerp = CGAL::ORIGIN + ((1 - u) * pos[v0] + u * pos[v1]);
+
+            const auto& p0 = domain.position(v0);
+            const auto& p1 = domain.position(v1);
+
+            if ((val0 <= iso_value) != (val1 <= iso_value)) {  // this edge is intersected by the isosurface
+                const FT u = (val0 - iso_value) / (val0 - val1);
+                const Point p_lerp = CGAL::ORIGIN + ((1 - u) * (p0 - CGAL::ORIGIN) + u * (p1 - CGAL::ORIGIN));
                 edge_intersections.push_back(p_lerp);
-                // const Vector n_lerp = (1 - u) * normals[v0] + u * normals[v1];
                 edge_intersection_normals.push_back(domain.gradient(p_lerp));
             }
         }
 
-        // MC Polygon Center of Mass
+        if (edge_intersections.empty()) {
-        if (false) {
+            return false;
-            Vector com_vec(0, 0, 0);
-
-            for (int i = 0; i < edge_intersections.size(); ++i) {
-                com_vec += edge_intersections[i] - CGAL::ORIGIN;
-            }
-
-            Point p = CGAL::ORIGIN + com_vec / edge_intersections.size();
-            point = p;
         }
 
         // SVD QEM
-        if (true) {
+        Eigen::Matrix3d A;
-            Eigen::Matrix3d A;
+        A.setZero();
-            A.setZero();
+        Eigen::Vector3d rhs;
-            Eigen::Vector3d b;
+        rhs.setZero();
-            b.setZero();
+        for (std::size_t i = 0; i < edge_intersections.size(); ++i) {
-            for (std::size_t i = 0; i < edge_intersections.size(); ++i) {
+            Eigen::Vector3d n_k = {edge_intersection_normals[i].x(), edge_intersection_normals[i].y(),
-                Eigen::Vector3d n_k = {edge_intersection_normals[i].x(), edge_intersection_normals[i].y(),
+                                   edge_intersection_normals[i].z()};
-                                       edge_intersection_normals[i].z()};
+            Eigen::Vector3d p_k = {edge_intersections[i].x(), edge_intersections[i].y(), edge_intersections[i].z()};
-                Eigen::Vector3d p_k = {edge_intersections[i].x(), edge_intersections[i].y(), edge_intersections[i].z()};
+            double d_k = n_k.transpose() * p_k;
-                double d_k = n_k.transpose() * p_k;
 
-                Eigen::Matrix3d A_k = n_k * n_k.transpose();
+            Eigen::Matrix3d A_k = n_k * n_k.transpose();
-                Eigen::Vector3d b_k = d_k * n_k;
+            Eigen::Vector3d b_k = d_k * n_k;
-                A += A_k;
+            A += A_k;
-                b += b_k;
+            rhs += b_k;
-            }
-
-            Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
-            // set threshold as in Peter Lindstrom's paper, "Out-of-Core
-            // Simplification of Large Polygonal Models"
-            svd.setThreshold(1e-3);
-
-            // Init x hat
-            Eigen::Vector3d x_hat;
-            x_hat << point.x(), point.y(), point.z();
-
-            // Lindstrom formula for QEM new position for singular matrices
-            Eigen::VectorXd v_svd = x_hat + svd.solve(b - A * x_hat);
-
-            point = Point(v_svd[0], v_svd[1], v_svd[2]);
         }
 
+        Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
+        // set threshold as in Peter Lindstrom's paper, "Out-of-Core
+        // Simplification of Large Polygonal Models"
+        svd.setThreshold(1e-3);
+
+        // Init x hat
+        Eigen::Vector3d x_hat;
+        x_hat << point.x(), point.y(), point.z();
+
+        // Lindstrom formula for QEM new position for singular matrices
+        Eigen::VectorXd v_svd = x_hat + svd.solve(rhs - A * x_hat);
+
+        point = Point(v_svd[0], v_svd[1], v_svd[2]);
+
         // bbox
         if (use_bbox) {
-            CGAL::Bbox_3 bbox = (CGAL::ORIGIN + pos[0]).bbox() + (CGAL::ORIGIN + pos[7]).bbox();
+            CGAL::Bbox_3 bbox = pos[0].bbox() + pos[7].bbox();  // TODO remove[0],[7]
 
             FT x = std::min<FT>(std::max<FT>(point.x(), bbox.xmin()), bbox.xmax());
             FT y = std::min<FT>(std::max<FT>(point.y(), bbox.ymin()), bbox.ymax());
@@ -155,123 +146,119 @@ public:
     }
 };
 
-class Voxel_center {
+/**
+ * \ingroup PkgIsosurfacing3Ref
+ *
+ * \brief Returns cell center.
+ */
+class Cell_center {
 public:
-    /// <summary>
+    /**
-    /// Compute vertex position for Dual Contouring
+     * \ingroup PkgIsosurfacing3Ref
-    /// </summary>
+     *
-    /// <typeparam name="Domain_"></typeparam>
+     * \brief Computes the vertex position for a point in Dual Contouring.
-    /// <param name="domain"></param>
+     *
-    /// <param name="iso_value"></param>
+     * \details
-    /// <param name="i"></param>
+     *
-    /// <param name="j"></param>
+     * \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
-    /// <param name="k"></param>
+     *
-    /// <returns> true, if there is a point in the cell</returns>
+     * \param domain the domain providing input data and its topology
+     * \param iso_value value of the isosurface
+     * \param cell the cell within the domain for which the vertex position ins computed
+     * \param point the point position of the vertex that belongs to that cell
+     *
+     * \return true, if the voxel intersects the isosurface
+     */
     template <class Domain_>
     bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh,
                   typename Domain_::Point& point) const {
         typedef typename Domain_::Point Point;
-        typedef typename Domain_::Vector_3 Vector;
+        typedef typename Domain_::Geom_traits::Vector_3 Vector;
+        typedef typename Domain_::FT FT;
 
-        namespace Tables = internal::Cube_table;
+        typename Domain_::Cell_vertices vertices = domain.cell_vertices(vh);
 
-        std::array<typename Domain_::FT, Tables::N_VERTICES> s = domain.voxel_values(vh);
+        std::vector<Point> pos(vertices.size());
+        std::transform(vertices.begin(), vertices.end(), pos.begin(),
+                       [&](const auto& v) { return domain.position(v); });
 
-        std::array<bool, Tables::N_VERTICES> b;
+        point = CGAL::centroid(pos.begin(), pos.end(), CGAL::Dimension_tag<0>());  // set point to cell center
-        std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; });
 
-        unsigned int cubeindex = 0;
+        bool allSmaller = true;
-        // set bit if corresponding corner is below iso
+        bool allGreater = true;
-        for (int i = 0; i < Tables::N_VERTICES; ++i) {
+        for (const auto& v : vertices) {
-            cubeindex |= b[i] << i;
+            const bool& b = domain.value(v) <= iso_value;
+            allSmaller = allSmaller && b;
+            allGreater = allGreater && !b;
         }
 
-        if (cubeindex == 0 || cubeindex == 255) {
+        if (allSmaller || allGreater) {
             return false;
         }
 
-        std::array<Point, Tables::N_VERTICES> p = domain.voxel_vertex_positions(vh);
-        std::array<Vector, Tables::N_VERTICES> pos;
-        std::transform(p.begin(), p.end(), pos.begin(), [](const auto& e) { return e - CGAL::ORIGIN; });
-
-        point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0]));  // set point to voxel center
-
         return true;
     }
 };
 
-class MC_polygon_center {
+/**
+ * \ingroup PkgIsosurfacing3Ref
+ *
+ * \brief Computes the centroid of all cell edge intersections with the isosurface.
+ */
+class Centroid_of_edge_intersections {
 public:
-    /// <summary>
+    /**
-    /// Compute vertex position for Dual Contouring
+     * \ingroup PkgIsosurfacing3Ref
-    /// </summary>
+     *
-    /// <typeparam name="Domain_"></typeparam>
+     * \brief Computes the vertex position for a point in Dual Contouring.
-    /// <param name="domain"></param>
+     *
-    /// <param name="iso_value"></param>
+     * \details
-    /// <param name="i"></param>
+     *
-    /// <param name="j"></param>
+     * \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`.
-    /// <param name="k"></param>
+     *
-    /// <returns> true, if there is a point in the cell</returns>
+     * \param domain the domain providing input data and its topology
+     * \param iso_value value of the isosurface
+     * \param cell the cell within the domain for which the vertex position ins computed
+     * \param point the point position of the vertex that belongs to that cell
+     *
+     * \return true, if the voxel intersects the isosurface
+     */
     template <class Domain_>
-    bool position(const Domain_& domain, const typename Domain_::FT iso_value, const typename Domain_::Cell_handle& vh,
+    bool position(const Domain_& domain, const typename Domain_::FT iso_value,
-                  typename Domain_::Point& point) const {
+                  const typename Domain_::Cell_handle& cell, typename Domain_::Point& point) const {
         typedef typename Domain_::Point Point;
-        typedef typename Domain_::Vector_3 Vector;
+        typedef typename Domain_::Geom_traits::Vector_3 Vector;
         typedef typename Domain_::FT FT;
 
-        namespace Tables = internal::Cube_table;
+        typename Domain_::Cell_vertices vertices = domain.cell_vertices(cell);
-
-        std::array<typename Domain_::FT, Tables::N_VERTICES> s = domain.voxel_values(vh);
-
-        std::array<bool, Tables::N_VERTICES> b;
-        std::transform(s.begin(), s.end(), b.begin(), [iso_value](const auto& e) { return e <= iso_value; });
-
-        unsigned int cubeindex = 0;
-        // set bit if corresponding corner is below iso
-        for (int i = 0; i < Tables::N_VERTICES; ++i) {
-            cubeindex |= b[i] << i;
-        }
-
-        if (cubeindex == 0 || cubeindex == 255) {
-            return false;
-        }
-
-        std::array<Point, Tables::N_VERTICES> p = domain.voxel_vertex_positions(vh);
-        std::array<Vector, Tables::N_VERTICES> pos;
-        std::transform(p.begin(), p.end(), pos.begin(), [](const auto& e) { return e - CGAL::ORIGIN; });
-
-        point = CGAL::ORIGIN + (pos[0] + 0.5 * (pos[7] - pos[0]));  // set point to voxel center
-
-        std::array<Vector, Tables::N_VERTICES> normals;
-        std::transform(vertices.begin(), vertices.end(), normals.begin(),
-                       [&](const auto& v) { return domain.gradient(domain.position(v)); });
-
 
         // compute edge intersections
         std::vector<Point> edge_intersections;
-        std::vector<Vector> edge_intersection_normals;
 
-        for (int i = 0; i < Tables::N_EDGES; ++i) {
+        for (const auto& edge : domain.cell_edges(cell)) {
-            const auto& v0 = Tables::edge_to_vertex[i][0];
+            const auto& edge_vertices = domain.edge_vertices(edge);
-            const auto& v1 = Tables::edge_to_vertex[i][1];
+            const auto& v0 = edge_vertices[0];
+            const auto& v1 = edge_vertices[1];
 
-            if (b[v0] != b[v1]) {  // e0
+            const auto& val0 = domain.value(v0);
-                const FT u = (s[v0] - iso_value) / (s[v0] - s[v1]);
+            const auto& val1 = domain.value(v1);
-                const Point p_lerp = CGAL::ORIGIN + ((1 - u) * pos[v0] + u * pos[v1]);
+
+            const auto& p0 = domain.position(v0);
+            const auto& p1 = domain.position(v1);
+
+            if ((val0 <= iso_value) != (val1 <= iso_value)) {  // this edge is intersected by the isosurface
+                const FT u = (val0 - iso_value) / (val0 - val1);
+                const Point p_lerp = CGAL::ORIGIN + ((1 - u) * (p0 - CGAL::ORIGIN) + u * (p1 - CGAL::ORIGIN));
                 edge_intersections.push_back(p_lerp);
-                const Vector n_lerp = (1 - u) * normals[v0] + u * normals[v1];
-                edge_intersection_normals.push_back(n_lerp);
             }
         }
 
-        // MC Polygon Center of Mass
+        if (edge_intersections.empty()) {
-        Vector com_vec(0, 0, 0);
+            return false;
-
-        for (int i = 0; i < edge_intersections.size(); ++i) {
-            com_vec += edge_intersections[i] - CGAL::ORIGIN;
         }
 
-        point = CGAL::ORIGIN + com_vec / edge_intersections.size();
+        point = CGAL::centroid(edge_intersections.begin(), edge_intersections.end(),
+                               CGAL::Dimension_tag<0>());  // set point to center of edge intersections
 
         return true;
     }
@@ -279,7 +266,7 @@ public:
 }  // namespace Positioning
 
 template <class Domain_, class Positioning_>
-class Dual_contouring_position_functor {
+class Dual_contouring_vertex_positioning {
 private:
     typedef Domain_ Domain;
     typedef Positioning_ Positioning;
@@ -289,7 +276,7 @@ private:
     typedef typename Domain::Cell_handle Cell_handle;
 
 public:
-    Dual_contouring_position_functor(const Domain& domain, FT iso_value, const Positioning& positioning)
+    Dual_contouring_vertex_positioning(const Domain& domain, FT iso_value, const Positioning& positioning)
         : domain(domain), iso_value(iso_value), positioning(positioning), points_counter(0) {}
 
     void operator()(const Cell_handle& v) {
@@ -316,7 +303,7 @@ public:
 };
 
 template <class Domain_>
-class Dual_contouring_quads_functor {
+class Dual_contouring_face_generation {
 private:
     typedef Domain_ Domain;
 
@@ -325,30 +312,30 @@ private:
     typedef typename Domain_::Cell_handle Cell_handle;
 
 public:
-    Dual_contouring_quads_functor(const Domain& domain, FT iso_value) : domain(domain), iso_value(iso_value) {}
+    Dual_contouring_face_generation(const Domain& domain, FT iso_value) : domain(domain), iso_value(iso_value) {}
 
     void operator()(const Edge_handle& e) {
-        // save all quads
+        // save all faces
         const auto& vertices = domain.edge_vertices(e);
         const FT s0 = domain.value(vertices[0]);
         const FT s1 = domain.value(vertices[1]);
 
         if (s0 <= iso_value && s1 > iso_value) {
-            const auto voxels = domain.cells_incident_to_edge(e);
+            const auto& voxels = domain.cells_incident_to_edge(e);
 
             std::lock_guard<std::mutex> lock(mutex);
-            quads[e] = {voxels[0], voxels[1], voxels[2], voxels[3]};
+            faces[e].insert(faces[e].begin(), voxels.begin(), voxels.end());
 
         } else if (s1 <= iso_value && s0 > iso_value) {
-            const auto voxels = domain.cells_incident_to_edge(e);
+            const auto& voxels = domain.cells_incident_to_edge(e);
 
             std::lock_guard<std::mutex> lock(mutex);
-            quads[e] = {voxels[0], voxels[3], voxels[2], voxels[1]};
+            faces[e].insert(faces[e].begin(), voxels.rbegin(), voxels.rend());
         }
     }
 
     // private:
-    std::map<Edge_handle, std::array<Cell_handle, 4>> quads;
+    std::map<Edge_handle, std::vector<Cell_handle>> faces;
 
     const Domain& domain;
     FT iso_value;

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h

@@ -43,7 +43,11 @@
 
 #include <CGAL/license/Isosurfacing_3.h>
 #include <CGAL/Isosurfacing_3/internal/Tables.h>
+#ifdef CGAL_LINKED_WITH_TBB
 #include <tbb/concurrent_vector.h>
+#else
+#include <vector>
+#endif
 
 #include <array>
 #include <atomic>
@@ -136,12 +140,7 @@ void mc_construct_triangles(const int i_case, const Vertices_& vertices, Triangl
         const int eg2 = Cube_table::triangle_cases[t_index + 2];
 
         // insert new triangle in list
-        std::array<Point, 3> points;
+        triangles.push_back({vertices[eg0], vertices[eg1], vertices[eg2]});
-        points[0] = vertices[eg0];
-        points[1] = vertices[eg1];
-        points[2] = vertices[eg2];
-
-        triangles.push_back(points);
     }
 }
 
@@ -154,11 +153,7 @@ void to_indexed_face_set(const TriangleList& triangle_list, PointRange& points,
         points.push_back(triangle[1]);
         points.push_back(triangle[2]);
 
-        polygons.push_back({});
+        polygons.push_back({id + 2, id + 1, id + 0});
-        auto& triangle = polygons.back();
-        triangle.push_back(id + 2);
-        triangle.push_back(id + 1);
-        triangle.push_back(id + 0);
     }
 }
 
@@ -170,7 +165,11 @@ private:
     typedef typename Domain::Point Point;
     typedef typename Domain::Cell_handle Cell_handle;
 
+#ifdef CGAL_LINKED_WITH_TBB
     typedef tbb::concurrent_vector<std::array<Point, 3>> Triangle_list;
+#else
+    typedef std::vector<std::array<Point, 3>> Triangle_list;
+#endif
 
 public:
     Marching_cubes_functor(const Domain& domain, const FT iso_value) : domain(domain), iso_value(iso_value) {}

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Tmc_internal.h

@@ -45,6 +45,7 @@
 #include <CGAL/Isosurfacing_3/internal/Marching_cubes_3_internal.h>
 #include <CGAL/Isosurfacing_3/internal/Tables.h>
 
+#include <cmath>
 #include <array>
 #include <atomic>
 #include <map>

Isosurfacing_3/package_info/Isosurfacing_3/dependencies

@@ -1,15 +1,18 @@
 Algebraic_foundations
-BGL
+CGAL_ImageIO
-Circulator
 Distance_2
-Hash_map
+Distance_3
 Installation
+Intersections_2
+Intersections_3
 Interval_support
+Isosurfacing_3
 Kernel_23
-Iso_surfacing_3
-Mesher_level
 Modular_arithmetic
 Number_types
+Orthtree
+Principal_component_analysis_LGPL
 Profiling_tools
+Property_map
 STL_Extension
 Stream_support