From d4528dd067a6a07db475ff45b7a15837699804a0 Mon Sep 17 00:00:00 2001 From: Julian Stahl Date: Sat, 3 Dec 2022 20:12:38 +0100 Subject: [PATCH] Improve reference manual --- .../Concepts/IsosurfacingDomain.h | 57 ++++++++-- .../Concepts/IsosurfacingDomainWithGradient.h | 7 +- .../doc/Isosurfacing_3/examples.txt | 2 +- .../include/CGAL/Cartesian_grid_3.h | 105 +++++++++++++++++- .../include/CGAL/Dual_contouring_3.h | 20 ++-- .../CGAL/Explicit_cartesian_grid_domain.h | 14 ++- .../CGAL/Implicit_cartesian_grid_domain.h | 26 ++++- .../include/CGAL/Marching_cubes_3.h | 24 ++-- 8 files changed, 212 insertions(+), 43 deletions(-) diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h index bfd3cf706b2..c4e2543a576 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h +++ b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h @@ -7,7 +7,6 @@ fulfilled by any class used as input data for any isosurfacing algorithms. \cgalHasModel `CGAL::Isosurfacing::Explicit_cartesian_grid_domain` \cgalHasModel `CGAL::Isosurfacing::Implicit_cartesian_grid_domain` -\cgalHasModel `CGAL::Isosurfacing::Implicit_octree_domain` */ @@ -48,21 +47,25 @@ public: /*! A container for the two vertices of an edge. + Must be a model of the concept `RandomAccessContainer` with size 2 whose value type is `Vertex_descriptor`. */ typedef unspecified_type Vertices_incident_to_edge; /*! A container for the cells incident to an edge. + Must be a model of the concept `RandomAccessContainer` whose value type is `Cell_descriptor`. */ typedef unspecified_type Cells_incident_to_edge; /*! A container for the vertices of a cell. + Must be a model of the concept `RandomAccessContainer` whose value type is `Vertex_descriptor`. */ typedef unspecified_type Cell_vertices; /*! A container for the edges of a cell. + Must be a model of the concept `RandomAccessContainer` whose value type is `Edge_descriptor`. */ typedef unspecified_type Cell_edges; @@ -74,67 +77,103 @@ public: /// @{ /*! - Returns the position of vertex v in 3D space + Get the position of a vertex in 3D space + + \param v the descriptor of the vertex + + \return the position of the vertex as a point */ Point position(const Vertex_descriptor& v) const; /*! - Returns the value of vertex v + Get the value of the function at a vertex + + \param v the descriptor of the vertex + + \return the value of the function */ FT value(const Vertex_descriptor& v) const; /*! - Returns the two vertices incident to edge e + Get the two vertices incident to an edge + + \param e the descriptor of the edge + + \return a collection of the two vertex descriptors */ Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const; /*! - Returns all voxels incident to edge e + Get all cells incident to an edge + + \param e the descriptor of the edge + + \return a collection of cell descriptors */ Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const; /*! - Returns all vertices of the cell c + Get all vertices of the a cell + + \param c the descriptor of the cell + + \return a collection of vertex descriptors */ Cell_vertices cell_vertices(const Cell_descriptor& c) const; /*! - Returns all edges of the cell c + Get all edges of the cell c + + \param c the descriptor of the cell + + \return a collection of edge descriptors */ Cell_edges cell_edges(const Cell_descriptor& c) const; /*! - Iterate sequentially over all vertices and call the functor f on each one + Iterate sequentially over all vertices and call a functor on each one + + \param f the functor called with every vertex. Must implement `void operator()(const Vertex_descriptor& vertex)`. */ template void iterate_vertices(Functor& f, Sequential_tag) const; /*! Iterate sequentially over all edges and call the functor f on each one + + \param f the functor called with every edge. Must implement `void operator()(const Edge_descriptor& edge)`. */ template void iterate_edges(Functor& f, Sequential_tag) const; /*! Iterate sequentially over all cells and call the functor f on each one + + \param f the functor called with every cell. Must implement `void operator()(const Cell_descriptor& cell)`. */ template void iterate_cells(Functor& f, Sequential_tag) const; /*! - (Optional) Iterate in parallel over all vertices and call the functor f on each one + (Optional) Iterate in parallel over all vertices and call a functor on each one + + \param f the functor called with every vertex. Must implement `void operator()(const Vertex_descriptor& vertex)`. */ template void iterate_vertices(Functor& f, Parallel_tag) const; /*! (Optional) Iterate in parallel over all edges and call the functor f on each one + + \param f the functor called with every edge. Must implement `void operator()(const Edge_descriptor& edge)`. */ template void iterate_edges(Functor& f, Parallel_tag) const; /*! (Optional) Iterate in parallel over all cells and call the functor f on each one + + \param f the functor called with every cell. Must implement `void operator()(const Cell_descriptor& cell)`. */ template void iterate_cells(Functor& f, Parallel_tag) const; diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomainWithGradient.h b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomainWithGradient.h index 171cc46e55d..2d3e105efd5 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomainWithGradient.h +++ b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomainWithGradient.h @@ -7,7 +7,6 @@ fulfilled by any class used as input data for some isosurfacing algorithms. \cgalHasModel `CGAL::Isosurfacing::Explicit_cartesian_grid_domain` \cgalHasModel `CGAL::Isosurfacing::Implicit_cartesian_grid_domain` -\cgalHasModel `CGAL::Isosurfacing::Implicit_octree_domain` */ @@ -28,7 +27,11 @@ public: /// @{ /*! - Returns the gradient at the position `p` + Get the gradient at a position + + \param p the point at which the gradient is evaluated + + \return the gradient vector */ Vector gradient(const Point& p) const; diff --git a/Isosurfacing_3/doc/Isosurfacing_3/examples.txt b/Isosurfacing_3/doc/Isosurfacing_3/examples.txt index b5b64b2a4e0..28df4609397 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/examples.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/examples.txt @@ -2,7 +2,7 @@ \example Isosurfacing_3/marching_cubes_implicit_sphere.cpp \example Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp \example Isosurfacing_3/marching_cubes_inrimage.cpp -\example Isosurfacing_3/marching_cubes_mesh_offset.cpp +\example Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp \example Isosurfacing_3/dual_contouring_cartesian_grid.cpp \example Isosurfacing_3/dual_contouring_mesh_offset.cpp \example Isosurfacing_3/dual_contouring_octree.cpp diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h index 84b42ecf069..2a012dd5f00 100644 --- a/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h +++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_3.h @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -23,6 +24,13 @@ namespace CGAL { +/** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Stores scalar values and gradients at points of a cartesian grid. + * + * \tparam GeomTraits the traits type + */ template class Cartesian_grid_3 { public: @@ -30,31 +38,87 @@ public: typedef typename Geom_traits::FT FT; typedef typename Geom_traits::Vector_3 Vector; + typedef internal::Grid_topology::Vertex_descriptor VertexDescriptor; + public: + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Create a grid with `xdim * ydim * zdim` grid points. + * The grid covers the space described by a bounding box. + * + * \param xdim the number of grid points in x direction + * \param ydim the number of grid points in y direction + * \param zdim the number of grid points in z direction + * \param bbox the bounding box + */ Cartesian_grid_3(const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, const Bbox_3& bbox) : sizes{xdim, ydim, zdim}, bbox(bbox) { + // pre-allocate memory values.resize(xdim * ydim * zdim); gradients.resize(xdim * ydim * zdim); + // calculate grid spacing const FT d_x = bbox.x_span() / (xdim - 1); const FT d_y = bbox.y_span() / (ydim - 1); const FT d_z = bbox.z_span() / (zdim - 1); spacing = Vector(d_x, d_y, d_z); } + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Create a grid from an `Image_3`. + * The dimensions and bounding box are read from the image. + * The values stored in the image must be of type `Geom_traits::FT` or implicitly convertible to it. + * + * \param image the image providing the data + */ Cartesian_grid_3(const Image_3& image) { from_image(image); } - FT operator()(const std::array& idx) const { - return values[linear_index(idx[0], idx[1], idx[2])]; + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Get the value stored in the grid at the grid point described by the vertex. + * + * \param v the vertex descriptor of the vertex + * + * \return the stored value + */ + FT operator()(const VertexDescriptor& v) const { + return values[linear_index(v[0], v[1], v[2])]; } + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Get the value stored in the grid at the given indices. + * + * \param x the index in x direction + * \param y the index in y direction + * \param z the index in z direction + * + * \return the stored value + */ FT value(const std::size_t x, const std::size_t y, const std::size_t z) const { return values[linear_index(x, y, z)]; } + /** + * \ingroup PkgIsosurfacing3Ref + * + * \brief Get a reference to the value stored in the grid at the given indices. + * Can be used to set values of the grid. + * + * \param x the index in x direction + * \param y the index in y direction + * \param z the index in z direction + * + * \return a reference to the stored value + */ FT& value(const std::size_t x, const std::size_t y, const std::size_t z) { return values[linear_index(x, y, z)]; } @@ -67,12 +131,29 @@ public: return gradients[linear_index(x, y, z)]; } + /** + * \ingroup PkgIsosurfacing3Ref + * + * \return the number of grid points in x direction + */ std::size_t xdim() const { return sizes[0]; } + + /** + * \ingroup PkgIsosurfacing3Ref + * + * \return the number of grid points in y direction + */ std::size_t ydim() const { return sizes[1]; } + + /** + * \ingroup PkgIsosurfacing3Ref + * + * \return the number of grid points in z direction + */ std::size_t zdim() const { return sizes[2]; } @@ -87,6 +168,7 @@ public: private: std::size_t linear_index(const std::size_t x, const std::size_t y, const std::size_t z) const { + // convert x, y, z into a linear index to access the vector return (z * ydim() + y) * xdim() + x; } @@ -105,20 +187,25 @@ private: template void Cartesian_grid_3::from_image(const Image_3& image) { + // compute bounding box const FT max_x = image.tx() + (image.xdim() - 1) * image.vx(); const FT max_y = image.ty() + (image.ydim() - 1) * image.vy(); const FT max_z = image.tz() + (image.zdim() - 1) * image.vz(); bbox = Bbox_3(image.tx(), image.ty(), image.tz(), max_x, max_y, max_z); + // get spacing spacing = Vector(image.vx(), image.vy(), image.vz()); + // get sizes sizes[0] = image.xdim(); sizes[1] = image.ydim(); sizes[2] = image.zdim(); + // pre-allocate values.resize(xdim() * ydim() * zdim()); gradients.resize(xdim() * ydim() * zdim()); + // copy values for (std::size_t x = 0; x < sizes[0]; x++) { for (std::size_t y = 0; y < sizes[1]; y++) { for (std::size_t z = 0; z < sizes[2]; z++) { @@ -131,23 +218,26 @@ void Cartesian_grid_3::from_image(const Image_3& image) { template Image_3 Cartesian_grid_3::to_image() const { - + // select the number type WORD_KIND wordkind; if (std::is_floating_point::value) wordkind = WK_FLOAT; else wordkind = WK_FIXED; + // select signed or unsigned SIGN sign; if (std::is_signed::value) sign = SGN_SIGNED; else sign = SGN_UNSIGNED; - const double vx = bbox.x_span() / (xdim() - 1); - const double vy = bbox.y_span() / (ydim() - 1); - const double vz = bbox.z_span() / (zdim() - 1); + // get spacing + const double vx = spacing()[0]; + const double vy = spacing()[1]; + const double vz = spacing()[2]; + // create image _image* im = _createImage(xdim(), ydim(), zdim(), 1, // vectorial dimension vx, vy, vz, // voxel size @@ -155,14 +245,17 @@ Image_3 Cartesian_grid_3::to_image() const { wordkind, // image word kind WK_FIXED, WK_FLOAT, WK_UNKNOWN sign); // image word sign + // error handling if (im == nullptr || im->data == nullptr) { throw std::bad_alloc(); // TODO: idk? } + // set min coordinates im->tx = bbox.xmin(); im->ty = bbox.ymin(); im->tz = bbox.zmin(); + // copy data FT* data = (FT*)im->data; for (std::size_t x = 0; x < xdim(); x++) { for (std::size_t y = 0; y < ydim(); y++) { diff --git a/Isosurfacing_3/include/CGAL/Dual_contouring_3.h b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h index ed59b038316..aaf53025bb0 100644 --- a/Isosurfacing_3/include/CGAL/Dual_contouring_3.h +++ b/Isosurfacing_3/include/CGAL/Dual_contouring_3.h @@ -15,7 +15,6 @@ #include -#include #include #include @@ -29,12 +28,12 @@ namespace Isosurfacing { * * \details * - * \tparam ConcurrencyTag determines if the algorithm is executed sequentially or in parallel. + * \tparam ConcurrencyTag determines if the algorithm is executed sequentially or in parallel. Default is sequential. * - * \tparam Domain_ must be a model of `IsosurfacingDomain`. + * \tparam Domain_ must be a model of `IsosurfacingDomainWithGradient`. * * \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. + * be constructed from the point type of the domain. * * \tparam PolygonRange a model of the concept * `RandomAccessContainer` and `BackInsertionSequence` whose value type is itself a model of the concepts @@ -45,37 +44,40 @@ namespace Isosurfacing { * * \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 points points of the polzgons in the created 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 + * \param positioning the functor dealing with vertex positioning inside a cell */ template > void dual_contouring(const Domain_& domain, const typename Domain_::FT iso_value, PointRange& points, PolygonRange& polygons, const Positioning& positioning = Positioning()) { - // static_assert(Domain_::CELL_TYPE & ANY_CELL); - + // create vertices in each relevant cell internal::Dual_contouring_vertex_positioning pos_func(domain, iso_value, positioning); domain.template iterate_cells(pos_func); + // connect vertices around an edge to form a face internal::Dual_contouring_face_generation face_generation(domain, iso_value); domain.template iterate_edges(face_generation); - // write points and faces in ranges + // copy vertices to point range 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; } + // copy faces to polygon range polygons.reserve(face_generation.faces.size()); for (const auto& q : face_generation.faces) { std::vector vertex_ids; for (const auto& v_id : q.second) { + // ignore voxels that are outside the valid region and are not stored in the map 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]); } } + // ignore degenerated faces if (vertex_ids.size() > 2) { polygons.push_back(vertex_ids); } diff --git a/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_domain.h index 4e7fee4f7f8..c0798ba14d5 100644 --- a/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_domain.h +++ b/Isosurfacing_3/include/CGAL/Explicit_cartesian_grid_domain.h @@ -15,7 +15,7 @@ #include #include -#include +#include #include #include #include @@ -23,6 +23,15 @@ namespace CGAL { namespace Isosurfacing { +/** + * \ingroup PkgIsosurfacing3Ref + * + * \brief A domain that respesents an explicitly stored cartesian grid. It is a model of the concept + * `IsosurfacingDomainWithGradient`. + * + * \tparam GeomTraits the traits type + * \tparam Gradient_ the type of the gradient functor. It must implement `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`. + */ template using Explicit_cartesian_grid_domain = Base_domain, Cartesian_grid_3, Gradient_>; @@ -38,6 +47,8 @@ using Explicit_cartesian_grid_domain = Base_domain> Explicit_cartesian_grid_domain create_explicit_cartesian_grid_domain( @@ -57,6 +68,7 @@ Explicit_cartesian_grid_domain create_explicit_cartesian_ const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin()); const typename GeomTraits::Vector_3 spacing = grid->get_spacing(); + // create copies as shared_ptr for safe memory management const Topology topo = std::make_shared(size_i, size_j, size_k); const Geometry geom = std::make_shared(offset, spacing); const Function func = grid; diff --git a/Isosurfacing_3/include/CGAL/Implicit_cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Implicit_cartesian_grid_domain.h index df7083e16d6..7b20ab1bd80 100644 --- a/Isosurfacing_3/include/CGAL/Implicit_cartesian_grid_domain.h +++ b/Isosurfacing_3/include/CGAL/Implicit_cartesian_grid_domain.h @@ -14,7 +14,7 @@ #include -#include +#include #include #include #include @@ -23,6 +23,17 @@ namespace CGAL { namespace Isosurfacing { +/** + * \ingroup PkgIsosurfacing3Ref + * + * \brief A domain that respesents a cartesian grid that discretizes an implicit function. It is a model of the concept + * `IsosurfacingDomainWithGradient`. + * + * \tparam GeomTraits the traits type + * \tparam PointFunction the type of the implicit function. It must implement `GeomTraits::FT operator()(const + * GeomTraits::Point& point) const`. \tparam Gradient_ the type of the gradient functor. It must implement + * `GeomTraits::Vector operator()(const GeomTraits::Point& point) const`. + */ template using Implicit_cartesian_grid_domain = Base_domain, @@ -32,17 +43,22 @@ using Implicit_cartesian_grid_domain = /** * \ingroup PkgIsosurfacing3Ref * - * \brief Creates a domain from a `Cartesian_grid_3` that can be used as input for isosurfacing algorithms. + * \brief Creates a domain from an implicit function that can be used as input for isosurfacing algorithms. * - * \details + * \details The implicit function will be evaluated on the grid points of the virtual grid + * defined by the bounding box and spacing. By not storing any function values implicitly + * less memory accesses are required in comparison to an `Explicit_cartesian_grid_domain`. * * \tparam GeomTraits the traits type - * \tparam PointFunction the type of the implicit function + * \tparam PointFunction the type of the implicit function. It must implement `GeomTraits::FT operator()(const + * GeomTraits::Point& point) const`. * * \param bbox a bounding box that specifies the size of the functions domain * \param spacing the distance between discretized points on the function * \param point_function the function with a point as argument * \param gradient a function that describes the gradient of the data + * + * \return a new `Implicit_cartesian_grid_domain` */ template > Implicit_cartesian_grid_domain create_implicit_cartesian_grid_domain( @@ -56,12 +72,14 @@ Implicit_cartesian_grid_domain create_impl typedef typename Domain::Gradient Gradient; typedef typename Function::element_type::Point_function Point_function; + // calculate grid dimensions const std::size_t size_i = std::ceil(bbox.x_span() / spacing.x()) + 1; const std::size_t size_j = std::ceil(bbox.y_span() / spacing.y()) + 1; const std::size_t size_k = std::ceil(bbox.z_span() / spacing.z()) + 1; const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin()); + // create copies as shared_ptr for safe memory management const Topology topo = std::make_shared(size_i, size_j, size_k); const Geometry geom = std::make_shared(offset, spacing); const Point_function point_func = std::make_shared(point_function); diff --git a/Isosurfacing_3/include/CGAL/Marching_cubes_3.h b/Isosurfacing_3/include/CGAL/Marching_cubes_3.h index cad6e8a91a5..6451a3e5808 100644 --- a/Isosurfacing_3/include/CGAL/Marching_cubes_3.h +++ b/Isosurfacing_3/include/CGAL/Marching_cubes_3.h @@ -14,7 +14,6 @@ #include -#include #include #include @@ -24,38 +23,41 @@ namespace Isosurfacing { /** * \ingroup PkgIsosurfacing3Ref * - * \brief Creates a polygon soup that represents an isosurface using the marching cubes algorithm. + * \brief Creates a triangular indexed face set that represents an isosurface using the marching cubes algorithm. * * \details * - * \tparam ConcurrencyTag determines if the algorithm is executed sequentially or in parallel. + * \tparam ConcurrencyTag determines if the algorithm is executed sequentially or in parallel. Default is sequential. * * \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 + * be constructed from the point type of the domain. + * + * \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 + * \param points points of the triangles in the created indexed face set + * \param triangles each element in the vector describes a triangle using the indices of the points in `points` * \param topologically_correct decides whether the topologically correct variant of Marching Cubes should be used */ template void marching_cubes(const Domain_& domain, const typename Domain_::FT iso_value, PointRange& points, - TriangleRange& polygons, bool topologically_correct = true) { - - // static_assert(Domain_::CELL_TYPE & CUBICAL_CELL); + TriangleRange& triangles, bool topologically_correct = true) { if (topologically_correct) { - internal::TMC_functor functor(domain, iso_value, points, polygons); + // run TMC and directly write the result to points and triangles + internal::TMC_functor functor(domain, iso_value, points, triangles); domain.template iterate_cells(functor); } else { + // run MC internal::Marching_cubes_functor functor(domain, iso_value); domain.template iterate_cells(functor); - internal::to_indexed_face_set(functor.get_triangles(), points, polygons); + // copy the result to points and triangles + internal::to_indexed_face_set(functor.get_triangles(), points, triangles); } }