From 28b865194f1d0aaffeb69eed441e968ac231abb9 Mon Sep 17 00:00:00 2001 From: Julian Stahl Date: Thu, 18 Aug 2022 15:04:51 +0200 Subject: [PATCH] Add first part of reference manual --- .../Concepts/IsosurfacingDomain.h | 142 ++++++++++++ .../Concepts/MarchingCubeOracle.h | 32 --- .../doc/Isosurfacing_3/Isosurfacing_3.txt | 4 +- .../doc/Isosurfacing_3/PackageDescription.txt | 9 +- .../marching_cubes_cartesian_grid_sphere.cpp | 12 +- .../marching_cubes_implicit_sphere.cpp | 10 +- .../include/CGAL/Cartesian_grid_domain.h | 214 ++++++++++++++++++ Isosurfacing_3/include/CGAL/Marching_cubes.h | 39 ++++ Isosurfacing_3/include/CGAL/marching_cube.h | 0 9 files changed, 406 insertions(+), 56 deletions(-) create mode 100644 Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h delete mode 100644 Isosurfacing_3/doc/Isosurfacing_3/Concepts/MarchingCubeOracle.h create mode 100644 Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h create mode 100644 Isosurfacing_3/include/CGAL/Marching_cubes.h delete mode 100644 Isosurfacing_3/include/CGAL/marching_cube.h diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h new file mode 100644 index 00000000000..270f4cdceb9 --- /dev/null +++ b/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 + void iterate_vertices(Functor& f, Sequential_tag) const; + + /*! + Iterate sequentially over all edges and call the functor f on each one + */ + template + void iterate_edges(Functor& f, Sequential_tag) const; + + /*! + Iterate sequentially over all cells and call the functor f on each one + */ + template + void iterate_cells(Functor& f, Sequential_tag) const; + + /*! + (Optional) Iterate in parallel over all vertices and call the functor f on each one + */ + template + void iterate_vertices(Functor& f, Parallel_tag) const; + + /*! + (Optional) Iterate in parallel over all edges and call the functor f on each one + */ + template + void iterate_edges(Functor& f, Parallel_tag) const; + + /*! + (Optional) Iterate in parallel over all cells and call the functor f on each one + */ + template + void iterate_cells(Functor& f, Parallel_tag) const; + + /// @} +}; diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/MarchingCubeOracle.h b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/MarchingCubeOracle.h deleted file mode 100644 index 6c9c2fe980e..00000000000 --- a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/MarchingCubeOracle.h +++ /dev/null @@ -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; - - -/// @} - -}; - diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt b/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt index 700fee52ede..10a675ef889 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt +++ b/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. diff --git a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt b/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt index 795fe3bc1dc..49b6268cf3b 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt +++ b/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()` */ \ No newline at end of file diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp index a2a1a79b88c..7efdde2fedc 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp @@ -1,17 +1,13 @@ #include #include -#include #include -#include #include typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::Vector_3 Vector_3; typedef typename Kernel::Point_3 Point_3; -typedef CGAL::Surface_mesh Mesh; - typedef std::vector Point_range; typedef std::vector> 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); } \ No newline at end of file diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp index 40f45b977c0..2041233e074 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -1,9 +1,7 @@ #include #include -#include #include -#include #include typedef CGAL::Simple_cartesian 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); } \ No newline at end of file diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h new file mode 100644 index 00000000000..e46325d0a15 --- /dev/null +++ b/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 + +#include "Cartesian_grid_3.h" +#include "Isosurfacing_3/internal/Tables.h" + +namespace CGAL { +namespace Isosurfacing { + +template +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 Vertex_handle; + typedef std::array Edge_handle; + typedef std::array Cell_handle; + + typedef std::array Edge_vertices; + typedef std::array Cells_incident_to_edge; + typedef std::array Cell_vertices; + typedef std::array Cell_edges; + +public: + Cartesian_grid_domain(const Cartesian_grid_3& 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 + 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 + 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& 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(0, size_x - 1), iterator); + } +#endif // CGAL_LINKED_WITH_TBB + + template + 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 + 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& 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(0, size_x - 1), iterator); + } +#endif // CGAL_LINKED_WITH_TBB + + template + 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 + 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& 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(0, size_x - 1), iterator); + } +#endif // CGAL_LINKED_WITH_TBB + +private: + const Cartesian_grid_3* grid; +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_CARTESIAN_GRID_DOMAIN_H diff --git a/Isosurfacing_3/include/CGAL/Marching_cubes.h b/Isosurfacing_3/include/CGAL/Marching_cubes.h new file mode 100644 index 00000000000..59ca4f0f646 --- /dev/null +++ b/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 +void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typename Domain_::FT iso_value, + PointRange& points, PolygonRange& polygons) { + + internal::Marching_cubes_RG2 functor(domain, iso_value, points, polygons); + domain.iterate_cells(functor, ConcurrencyTag()); +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_MARCHING_CUBES_H diff --git a/Isosurfacing_3/include/CGAL/marching_cube.h b/Isosurfacing_3/include/CGAL/marching_cube.h deleted file mode 100644 index e69de29bb2d..00000000000