Add first part of reference manual

Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h

@@ -0,0 +1,142 @@
+/*!
+\ingroup PkgIsosurfacing3Concepts
+\cgalConcept
+
+The concept `IsosurfacingDomain` describes the set of requirements to be
+fulfilled by any class used as input data for any isosurfacing algorithms.
+
+\cgalHasModel `CGAL::Isosurfacing::Cartesian_grid_domain`
+\cgalHasModel `CGAL::Isosurfacing::Implicit_domain`
+
+*/
+
+class IsosurfacingDomain {
+public:
+    /// \name Types
+    /// @{
+
+    /*!
+    Traits type model of \cgal %Kernel
+    */
+    typedef unspecified_type Traits;
+
+    /*!
+    The scalar type.
+    */
+    typedef unspecified_type FT;
+
+    /*!
+    The point type.
+    */
+    typedef unspecified_type Point;
+
+    /*!
+    A handle to identify a vertex.
+    */
+    typedef unspecified_type Vertex_handle;
+
+    /*!
+    A handle to identify an edge.
+    */
+    typedef unspecified_type Edge_handle;
+
+    /*!
+    A handle to identify a cell.
+    */
+    typedef unspecified_type Cell_handle;
+
+    /*!
+    A container for the two vertices of an edge.
+    */
+    typedef unspecified_type Edge_vertices;
+
+    /*!
+    A container for the cells incident to an edge.
+    */
+    typedef unspecified_type Cells_incident_to_edge;
+
+    /*!
+    A container for the vertices of a cell.
+    */
+    typedef unspecified_type Cell_vertices;
+
+    /*!
+    A container for the edges of a cell.
+    */
+    typedef unspecified_type Cell_edges;
+
+
+    /// @}
+
+    /// \name Operations
+    /// The following member functions must exist.
+    /// @{
+
+    /*!
+    Returns the position of vertex v in 3D space
+    */
+    Point position(const Vertex_handle& v) const;
+
+    /*!
+    Returns the stored value of vertex v
+    */
+    FT value(const Vertex_handle& v) const;
+
+    /*!
+    Returns the two vertices incident to edge e
+    */
+    Edge_vertices edge_vertices(const Edge_handle& e) const;
+
+    /*!
+    Returns all voxels incident to edge e
+    */
+    Cells_incident_to_edge cells_incident_to_edge(const Edge_handle& e) const;
+
+    /*!
+    Returns all vertices of the cell c
+    */
+    Cell_vertices cell_vertices(const Cell_handle& c) const;
+
+    /*!
+    Returns all edges of the cell c
+    */
+    Cell_edges cell_edges(const Cell_handle& c) const;
+
+    /*!
+    Iterate sequentially over all vertices and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_vertices(Functor& f, Sequential_tag) const;
+
+    /*!
+    Iterate sequentially over all edges and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_edges(Functor& f, Sequential_tag) const;
+
+    /*!
+    Iterate sequentially over all cells and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_cells(Functor& f, Sequential_tag) const;
+
+    /*!
+    (Optional) Iterate in parallel over all vertices and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_vertices(Functor& f, Parallel_tag) const;
+
+    /*!
+    (Optional) Iterate in parallel over all edges and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_edges(Functor& f, Parallel_tag) const;
+
+    /*!
+    (Optional) Iterate in parallel over all cells and call the functor f on each one
+    */
+    template <typename Functor>
+    void iterate_cells(Functor& f, Parallel_tag) const;
+
+    /// @}
+};

Isosurfacing_3/doc/Isosurfacing_3/Concepts/MarchingCubeOracle.h

@@ -1,32 +0,0 @@
-/*!
-\ingroup PkgIsosurfacing3Concepts
-\cgalConcept
-
-The concept `MarchingCubeOracle` describes the set of requirements to be
-fulfilled by any class used by the marching cubes algorithms.>`.
-
-\cgalHasModel `CGAL::Marching_cubes_implicit_3`
-
-*/
-
-class MarchingCubeOracle {
-public:
-
-/// \name Types
-/// @{
-
-/*!
-The scalar type.
-*/
-typedef unspecified_type FT;
-
-/*!
-Traits type model of \cgal %Kernel
-*/
-typedef unspecified_type Traits;
-
-
-/// @}
-
-};
-

Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt

@@ -18,7 +18,7 @@ The representation that is used to store the isosurface is a triangular or quadr
 
 \section secmyalgorithms Algorithms
 
-There are multiple algorithms to extract an isosurface from a uniform grid.
+There are multiple algorithms to extract an isosurface from a uniform grid or octree.
 This package contains marching cubes, ...
 
 \section secmyinterface Interface
@@ -31,7 +31,7 @@ void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typena
                                              PointRange& points, PolygonRange& polygons);
 \endcode
 
-The `domain` is an object that provides functions to access the input data in the form of a uniform grid.
+The `domain` is an object that provides functions to access the input data and its topology.
 There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm.
 Examples are `Cartesian_grid_domain`, which takes a `Cartesian_grid_3` as parameter, 
 and `Implicit_domain`, which represents the grid as an implicit function.

Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt

@@ -3,7 +3,7 @@
 /// \ingroup PkgIsosurfacing3Ref
 
 /*!
-\addtogroup PkgIsosurfacingRef
+\addtogroup PkgIsosurfacing3Ref
 \cgalPkgDescriptionBegin{3D Isosurfacing,PkgIsosurfacing3}
 \cgalPkgPicture{isosurfacing3_ico.png}
 \cgalPkgSummaryBegin
@@ -24,14 +24,15 @@ Marching Cubes Algorithm .....
 \cgalClassifedRefPages
 
 \cgalCRPSection{Concepts}
-- `MarchingCubeOracle`
+- `IsosurfacingDomain`
 
 \cgalCRPSection{Classes}
 
-- `CGAL::Marching_cubes_implicit_3`
+- `CGAL::Isosurfacing::Cartesian_grid_domain`
+- `CGAL::Isosurfacing::Implicit_domain`
 
 \cgalCRPSection{Free Functions}
 
-- `CGAL::make_triangle_mesh_using_marching_cubes()`
+- `CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes()`
 
 */

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp

@@ -1,17 +1,13 @@
 
 #include <CGAL/Isosurfacing_3/Cartesian_grid_domain.h>
 #include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
-#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
 #include <CGAL/Simple_cartesian.h>
-#include <CGAL/Surface_mesh.h>
 #include <CGAL/boost/graph/IO/OFF.h>
 
 typedef CGAL::Simple_cartesian<double> Kernel;
 typedef typename Kernel::Vector_3 Vector_3;
 typedef typename Kernel::Point_3 Point_3;
 
-typedef CGAL::Surface_mesh<Point_3> Mesh;
-
 typedef std::vector<Point_3> Point_range;
 typedef std::vector<std::vector<std::size_t>> Polygon_range;
 
@@ -40,10 +36,6 @@ int main() {
     // execute marching cubes with an isovalue of 0.8
     CGAL::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
 
-    // convert the polygon soup to a surface mesh
-    Mesh mesh;
-    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
-
-    // save the mesh in the OFF format
-    CGAL::IO::write_OFF("result.off", mesh);
+    // save the polygon soup in the OFF format
+    CGAL::IO::write_OFF("result.off", points, polygons);
 }

Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp

@@ -1,9 +1,7 @@
 
 #include <CGAL/Isosurfacing_3/Implicit_domain.h>
 #include <CGAL/Isosurfacing_3/Marching_cubes_3.h>
-#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
 #include <CGAL/Simple_cartesian.h>
-#include <CGAL/Surface_mesh.h>
 #include <CGAL/boost/graph/IO/OFF.h>
 
 typedef CGAL::Simple_cartesian<double> Kernel;
@@ -32,10 +30,6 @@ int main() {
     // execute marching cubes with an isovalue of 0.8
     CGAL::make_triangle_mesh_using_marching_cubes(domain, 0.8f, points, polygons);
 
-    // convert the polygon soup to a surface mesh
-    Mesh mesh;
-    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh);
-
-    // save the mesh in the OFF format
-    CGAL::IO::write_OFF("result.off", mesh);
+    // save the polygon soup in the OFF format
+    CGAL::IO::write_OFF("result.off", points, polygons);
 }

Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h

@@ -0,0 +1,214 @@
+#ifndef CGAL_CARTESIAN_GRID_DOMAIN_H
+#define CGAL_CARTESIAN_GRID_DOMAIN_H
+
+#include <tbb/parallel_for.h>
+
+#include "Cartesian_grid_3.h"
+#include "Isosurfacing_3/internal/Tables.h"
+
+namespace CGAL {
+namespace Isosurfacing {
+
+template <class GeomTraits>
+class Cartesian_grid_domain {
+public:
+    typedef GeomTraits Geom_traits;
+    typedef typename Geom_traits::FT FT;
+    typedef typename Geom_traits::Point_3 Point;
+    typedef typename Geom_traits::Vector_3 Vector;
+
+    typedef std::array<std::size_t, 3> Vertex_handle;
+    typedef std::array<std::size_t, 4> Edge_handle;
+    typedef std::array<std::size_t, 3> Cell_handle;
+
+    typedef std::array<Vertex_handle, 2> Edge_vertices;
+    typedef std::array<Cell_handle, 4> Cells_incident_to_edge;
+    typedef std::array<Vertex_handle, internal::Cube_table::N_VERTICES> Cell_vertices;
+    typedef std::array<Edge_handle, internal::Cube_table::N_EDGES> Cell_edges;
+
+public:
+    Cartesian_grid_domain(const Cartesian_grid_3<Geom_traits>& grid) : grid(&grid) {}
+
+    Point position(const Vertex_handle& v) const {
+        const FT vx = grid->voxel_x();
+        const FT vy = grid->voxel_y();
+        const FT vz = grid->voxel_z();
+
+        return Point(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z());
+    }
+
+    Vector gradient(const Vertex_handle& v) const {
+        const FT vx = grid->voxel_x();
+        const FT vy = grid->voxel_y();
+        const FT vz = grid->voxel_z();
+
+        Vector g(v[0] * vx + grid->offset_x(), v[1] * vy + grid->offset_y(), v[2] * vz + grid->offset_z());
+        return g / std::sqrt(g.squared_length());
+    }
+
+    FT value(const Vertex_handle& v) const {
+        return grid->value(v[0], v[1], v[2]);
+    }
+
+    Edge_vertices edge_vertices(const Edge_handle& e) const {
+        Edge_vertices ev;
+        ev[0] = {e[0], e[1], e[2]};
+        ev[1] = {e[0], e[1], e[2]};
+        ev[1][e[3]] += 1;
+        return ev;
+    }
+
+    Cells_incident_to_edge cells_incident_to_edge(const Edge_handle& e) const {
+        const int local = internal::Cube_table::edge_store_index[e[3]];
+        auto neighbors = internal::Cube_table::edge_to_voxel_neighbor[local];
+
+        Cells_incident_to_edge cite;
+        for (std::size_t i = 0; i < cite.size(); i++) {
+            for (std::size_t j = 0; j < cite[i].size(); j++) {
+                cite[i][j] = e[j] + neighbors[i][j];
+            }
+        }
+        return cite;
+    }
+
+    Cell_vertices cell_vertices(const Cell_handle& v) const {
+        Cell_vertices cv;
+        for (std::size_t i = 0; i < cv.size(); i++) {
+            for (std::size_t j = 0; j < v.size(); j++) {
+                cv[i][j] = v[j] + internal::Cube_table::local_vertex_position[i][j];
+            }
+        }
+        return cv;
+    }
+
+    Cell_edges cell_edges(const Cell_handle& v) const {
+        Cell_edges ce;
+        for (std::size_t i = 0; i < ce.size(); i++) {
+            for (std::size_t j = 0; j < v.size(); j++) {
+                ce[i][j] = v[j] + internal::Cube_table::global_edge_id[i][j];
+            }
+            ce[i][3] = internal::Cube_table::global_edge_id[i][3];
+        }
+        return ce;
+    }
+
+    template <typename Functor>
+    void iterate_vertices(Functor& f, Sequential_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        for (std::size_t x = 0; x < size_x - 1; x++) {
+            for (std::size_t y = 0; y < size_y - 1; y++) {
+                for (std::size_t z = 0; z < size_z - 1; z++) {
+                    f({x, y, z});
+                }
+            }
+        }
+    }
+
+#ifdef CGAL_LINKED_WITH_TBB
+    template <typename Functor>
+    void iterate_vertices(Functor& f, Parallel_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        auto iterator = [=, &f](const tbb::blocked_range<std::size_t>& r) {
+            for (std::size_t x = r.begin(); x != r.end(); x++) {
+                for (std::size_t y = 0; y < size_y - 1; y++) {
+                    for (std::size_t z = 0; z < size_z - 1; z++) {
+                        f({x, y, z});
+                    }
+                }
+            }
+        };
+
+        tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_x - 1), iterator);
+    }
+#endif  // CGAL_LINKED_WITH_TBB
+
+    template <typename Functor>
+    void iterate_edges(Functor& f, Sequential_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        for (std::size_t x = 0; x < size_x - 1; x++) {
+            for (std::size_t y = 0; y < size_y - 1; y++) {
+                for (std::size_t z = 0; z < size_z - 1; z++) {
+                    f({x, y, z, 0});
+                    f({x, y, z, 1});
+                    f({x, y, z, 2});
+                }
+            }
+        }
+    }
+
+#ifdef CGAL_LINKED_WITH_TBB
+    template <typename Functor>
+    void iterate_edges(Functor& f, Parallel_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        auto iterator = [=, &f](const tbb::blocked_range<std::size_t>& r) {
+            for (std::size_t x = r.begin(); x != r.end(); x++) {
+                for (std::size_t y = 0; y < size_y - 1; y++) {
+                    for (std::size_t z = 0; z < size_z - 1; z++) {
+                        f({x, y, z, 0});
+                        f({x, y, z, 1});
+                        f({x, y, z, 2});
+                    }
+                }
+            }
+        };
+
+        tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_x - 1), iterator);
+    }
+#endif  // CGAL_LINKED_WITH_TBB
+
+    template <typename Functor>
+    void iterate_cells(Functor& f, Sequential_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        for (std::size_t x = 0; x < size_x - 1; x++) {
+            for (std::size_t y = 0; y < size_y - 1; y++) {
+                for (std::size_t z = 0; z < size_z - 1; z++) {
+                    f({x, y, z});
+                }
+            }
+        }
+    }
+
+#ifdef CGAL_LINKED_WITH_TBB
+    template <typename Functor>
+    void iterate_cells(Functor& f, Parallel_tag) const {
+        const std::size_t size_x = grid->xdim();
+        const std::size_t size_y = grid->ydim();
+        const std::size_t size_z = grid->zdim();
+
+        auto iterator = [=, &f](const tbb::blocked_range<std::size_t>& r) {
+            for (std::size_t x = r.begin(); x != r.end(); x++) {
+                for (std::size_t y = 0; y < size_y - 1; y++) {
+                    for (std::size_t z = 0; z < size_z - 1; z++) {
+                        f({x, y, z});
+                    }
+                }
+            }
+        };
+
+        tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_x - 1), iterator);
+    }
+#endif  // CGAL_LINKED_WITH_TBB
+
+private:
+    const Cartesian_grid_3<Geom_traits>* grid;
+};
+
+}  // namespace Isosurfacing
+}  // namespace CGAL
+
+#endif  // CGAL_CARTESIAN_GRID_DOMAIN_H

Isosurfacing_3/include/CGAL/Marching_cubes.h

@@ -0,0 +1,39 @@
+#ifndef CGAL_MARCHING_CUBES_H
+#define CGAL_MARCHING_CUBES_H
+
+#include "Isosurfacing_3/internal/Marching_cubes_3_internal.h"
+
+namespace CGAL {
+namespace Isosurfacing {
+
+/**
+ * \ingroup PkgIsosurfacing3Ref
+ *
+ * \brief Creates a polygon soup that represents an isosurface using the marching cubes 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`.
+ *
+ * \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 soup
+ * \param polygons each element in the vector describes a polygon using the indices of the points in points
+ */
+template <typename ConcurrencyTag = Sequential_tag, class Domain_, class PointRange, class PolygonRange>
+void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typename Domain_::FT iso_value,
+                                             PointRange& points, PolygonRange& polygons) {
+
+    internal::Marching_cubes_RG2<Domain_, PointRange, PolygonRange> functor(domain, iso_value, points, polygons);
+    domain.iterate_cells(functor, ConcurrencyTag());
+}
+
+}  // namespace Isosurfacing
+}  // namespace CGAL
+
+#endif  // CGAL_MARCHING_CUBES_H