From 8e0140e6412a45b8063332bc66e2ad6e82e8158a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 15 Feb 2024 10:34:43 +0100 Subject: [PATCH] New version of Isosurfacing_3: - Completely restructure the domain classes to separate what is spatial partitioning, and what is value/gradient field definition. - Improve DC edge-isosurfacing intersection, factorize code - Refactor DC implementation to make it easier to use new (better) oracles - Add concepts for these oracles, and document their models --- .../CGAL/Isosurfacing_3/Cartesian_grid_3.h | 461 ++++----- .../Isosurfacing_3/Dual_contouring_domain_3.h | 110 ++ .../Explicit_Cartesian_grid_domain_3.h | 107 -- .../Explicit_Cartesian_grid_gradient_3.h | 138 --- .../Finite_difference_gradient_3.h | 24 +- .../CGAL/Isosurfacing_3/Gradient_function_3.h | 90 ++ .../include/CGAL/Isosurfacing_3/IO/Image_3.h | 170 ++++ .../include/CGAL/Isosurfacing_3/IO/OBJ.h | 100 ++ .../Implicit_Cartesian_grid_domain_3.h | 172 ---- .../Isosurfacing_3/Implicit_octree_domain.h | 108 -- .../Interpolated_discrete_gradients_3.h | 124 +++ .../Interpolated_discrete_values_3.h | 124 +++ .../Isosurfacing_3/Marching_cubes_domain_3.h | 102 ++ .../CGAL/Isosurfacing_3/Value_function_3.h | 91 ++ .../CGAL/Isosurfacing_3/Zero_gradient.h | 46 - .../CGAL/Isosurfacing_3/dual_contouring_3.h | 90 +- .../edge_intersection_oracles_3.h | 212 ++++ .../CGAL/Isosurfacing_3/internal/Cell_type.h | 18 +- .../Explicit_Cartesian_grid_function.h | 52 - .../Explicit_Cartesian_grid_geometry_3.h | 47 - .../Implicit_Cartesian_grid_geometry_3.h | 67 -- .../Implicit_function_with_geometry.h | 53 - .../internal/Isosurfacing_domain_3.h | 135 +-- .../Isosurfacing_3/internal/Octree_geometry.h | 2 +- .../Isosurfacing_3/internal/Octree_topology.h | 2 +- .../Isosurfacing_3/internal/Octree_wrapper.h | 4 +- .../internal/dual_contouring_functors.h | 960 +++++++++++------- .../internal/implicit_shapes_helper.h | 182 ++++ .../internal/marching_cubes_functors.h | 38 +- .../internal/partition_traits.h | 26 + ..._3.h => partition_traits_Cartesian_grid.h} | 132 +-- .../internal/partition_traits_Octree.h | 27 + .../CGAL/Isosurfacing_3/internal/tables.h | 2 - ...ogically_correct_marching_cubes_functors.h | 132 ++- .../Isosurfacing_3/interpolation_schemes_3.h | 236 +++++ .../CGAL/Isosurfacing_3/marching_cubes_3.h | 15 +- 36 files changed, 2685 insertions(+), 1714 deletions(-) create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Dual_contouring_domain_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Gradient_function_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/OBJ.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Marching_cubes_domain_3.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Value_function_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/Zero_gradient.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/edge_intersection_oracles_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h delete mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits.h rename Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/{Grid_topology_3.h => partition_traits_Cartesian_grid.h} (60%) create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h create mode 100644 Isosurfacing_3/include/CGAL/Isosurfacing_3/interpolation_schemes_3.h diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h index a8b7c155745..ea834c18971 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Cartesian_grid_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,36 +8,71 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Julian Stahl +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_CARTESIAN_GRID_3_H #define CGAL_ISOSURFACING_3_CARTESIAN_GRID_3_H #include -#include +#include +#include #include #include #include #include -#include #include -#include #include -#include namespace CGAL { namespace Isosurfacing { /** - * \ingroup IS_Domains_grp + * \ingroup IS_Partitions_helpers_grp * - * \brief stores scalar values and gradients at the vertices of a %Cartesian grid. + * A policy to choose whether grid vertex positions should be cached, or recomputed at each access. + * + * \tparam Tag a tag that is either `Tag_true` (positions are cached) or `Tag_false` (positions are not cached). + */ +template +struct Grid_vertex_memory_policy : public Tag { }; + +/** + * \ingroup IS_Partitions_helpers_grp + * + * A convenience alias for the policy that caches grid vertex positions. + */ +using Cache_positions = Grid_vertex_memory_policy; + +/** + * \ingroup IS_Partitions_helpers_grp + * + * A convenience alias for the policy that does not cache grid vertex positions. + */ +using Do_not_cache_positions = Grid_vertex_memory_policy; + +/** + * \ingroup IS_Partitions_grp + * + * \cgalModels{Partition_3} + * + * \brief The class `Cartesian_grid_3` represents a 3D %Cartesian grid, that is the partition of + * an iso-cuboid into identical iso-cuboidal cells. + * + * The class `Cartesian_grid_3` is one of the possible space partitioning data structures + * that can be used along with value and gradient fields to make up a domain. * * \tparam GeomTraits must be a model of `IsosurfacingTraits_3`. + * \tparam MemoryPolicy whether the geometric positions of the grid vertices are stored or not. + * Possible values are `CGAL::Isosurfacing::Cache_positions` and `CGAL::Isosurfacing::Do_not_cache_positions`. + * + * \sa `CGAL::Isosurfacing::Marching_cubes_domain_3()` + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3()` */ -template +template // @todo actually implement it class Cartesian_grid_3 { public: @@ -49,51 +84,39 @@ public: private: Iso_cuboid_3 m_bbox; - std::array m_spacing; std::array m_sizes; - - std::vector m_values; - std::vector m_gradients; + std::array m_spacing; Geom_traits m_gt; public: /** - * \brief creates a grid with `xdim * ydim * zdim` grid vertices. + * \brief creates a %Cartesian grid with `xdim * ydim * zdim` grid vertices. * * The grid covers the space described by a bounding box. * + * \param bbox the bounding box of the grid * \param xdim the number of grid vertices in the `x` direction * \param ydim the number of grid vertices in the `y` direction * \param zdim the number of grid vertices in the `z` direction - * \param bbox the bounding box of the grid * \param gt the geometric traits * * \pre `xdim`, `ydim`, and `zdim` are (strictly) positive. */ - Cartesian_grid_3(const std::size_t xdim, + Cartesian_grid_3(const Iso_cuboid_3& bbox, + const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, - const Iso_cuboid_3& bbox, const Geom_traits& gt = Geom_traits()) - : m_sizes{xdim, ydim, zdim}, - m_bbox{bbox}, + : m_bbox{bbox}, + m_sizes{xdim, ydim, zdim}, m_gt{gt} { - CGAL_precondition(xdim > 0); - CGAL_precondition(ydim > 0); - CGAL_precondition(zdim > 0); - auto x_coord = m_gt.compute_x_3_object(); auto y_coord = m_gt.compute_y_3_object(); auto z_coord = m_gt.compute_z_3_object(); auto vertex = m_gt.construct_vertex_3_object(); - // pre-allocate memory - const std::size_t nv = xdim * ydim * zdim; - m_values.resize(nv); - m_gradients.resize(nv); - // calculate grid spacing const Point_3& min_p = vertex(bbox, 0); const Point_3& max_p = vertex(bbox, 7); @@ -107,58 +130,82 @@ public: m_spacing = make_array(d_x, d_y, d_z); } - Cartesian_grid_3(const std::size_t xdim, + /** + * \brief creates a %Cartesian grid with `xdim * ydim * zdim` grid vertices. + * + * The grid covers the space described by a bounding box, itself described through two diagonal corners. + * + * \param p the lowest corner of the bounding box of the grid + * \param q the upper corner of the bounding box of the grid + * \param xdim the number of grid vertices in the `x` direction + * \param ydim the number of grid vertices in the `y` direction + * \param zdim the number of grid vertices in the `z` direction + * \param gt the geometric traits + * + * \pre `p` is lexicographically (strictly) smaller than `q` + * \pre `xdim`, `ydim`, and `zdim` are (strictly) positive. + */ + Cartesian_grid_3(const Point_3& p, const Point_3& q, + const std::size_t xdim, const std::size_t ydim, const std::size_t zdim, - const Point_3& p, const Point_3& q, const Geom_traits& gt = Geom_traits()) - : Cartesian_grid_3(xdim, ydim, zdim, Iso_cuboid_3(p, q), gt) + : Cartesian_grid_3{Iso_cuboid_3{p, q}, xdim, ydim, zdim, gt} { } /** - * \brief creates a grid from a `CGAL::Image_3`. + * \brief creates a %Cartesian grid using a prescribed grid step. * - * 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. + * The grid covers the space described by a bounding box. * - * \param image the image providing the data + * \param bbox the bounding box of the grid + * \param spacing the dimension of the paving cell, in the `x`, `y`, and `z` directions, respectively. + * \param gt the geometric traits + * + * \pre the diagonal of `bbox` has length a multiple of `spacing` */ - Cartesian_grid_3(const Image_3& image) + Cartesian_grid_3(const Iso_cuboid_3& bbox, + const std::array& spacing, + const Geom_traits& gt = Geom_traits()) + : m_bbox{bbox}, + m_spacing{spacing}, + m_gt{gt} { - auto point = m_gt.construct_point_3_object(); - auto iso_cuboid = m_gt.construct_iso_cuboid_3_object(); + auto x_coord = gt.compute_x_3_object(); + auto y_coord = gt.compute_y_3_object(); + auto z_coord = gt.compute_z_3_object(); + auto vertex = gt.construct_vertex_3_object(); + auto vector = gt.construct_vector_3_object(); - // 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(); - m_bbox = iso_cuboid(point(image.tx(), image.ty(), image.tz()), - point(max_x, max_y, max_z)); + const Point_3& min_p = vertex(bbox, 0); + const Point_3& max_p = vertex(bbox, 7); + const FT x_span = x_coord(max_p) - x_coord(min_p); + const FT y_span = y_coord(max_p) - y_coord(min_p); + const FT z_span = z_coord(max_p) - z_coord(min_p); - // get spacing - m_spacing = make_array(image.vx(), image.vy(), image.vz()); - - // get sizes - m_sizes[0] = image.xdim(); - m_sizes[1] = image.ydim(); - m_sizes[2] = image.zdim(); - - // pre-allocate - const std::size_t nv = m_sizes[0] * m_sizes[1] * m_sizes[2]; - m_values.resize(nv); - m_gradients.resize(nv); - - // copy values - for(std::size_t x=0; x& spacing, + const Geom_traits& gt = Geom_traits()) + : Cartesian_grid_3{Iso_cuboid_3{p, q}, spacing, gt} + { } public: /** @@ -169,6 +216,17 @@ public: return m_gt; } + /** + * \return the bounding box of the %Cartesian grid + */ + const Iso_cuboid_3& bbox() const { return m_bbox; } + + /** + * \return the spacing of the %Cartesian grid, that is a vector whose coordinates are + * the grid steps in the `x`, `y`, and `z` directions, respectively + */ + const std::array& spacing() const { return m_spacing; } + /** * \return the number of grid vertices in the `x` direction */ @@ -184,251 +242,84 @@ public: */ std::size_t zdim() const { return m_sizes[2]; } +public: /** - * \return the bounding box of the %Cartesian grid. - */ - const Iso_cuboid_3& bbox() const { return m_bbox; } + * \brief gets the canonical index of a grid cell given its indices. + */ + std::size_t linear_index(const std::size_t i, + const std::size_t j, + const std::size_t k) const + { + CGAL_precondition(i < m_sizes[0] && j < m_sizes[1] && k < m_sizes[2]); - /** - * \return the spacing of the %Cartesian grid, that is a vector whose coordinates are - * the grid steps in the `x`, `y`, and `z` directions, respectively - */ - const std::array& spacing() const { return m_spacing; } + // convert (i, j, k) into a linear index, e.g. to access the scalar values / gradient vectors + return (k * m_sizes[1] + j) * m_sizes[0] + i; + } public: /** - * \brief gets the geometric position of the grid vertex described by a set of indices. + * \brief gets the index of the grid cell that contains a given point. * - * Positions are not stored but calculated from an offset and grid spacing. + * \param p the point to be located + * + * \return the index of the grid cell that contains `p` + * + * \pre `p` is inside the bounding box of the grid. + */ + std::array index(const Point_3& p) const + { + auto x_coord = m_gt.compute_x_3_object(); + auto y_coord = m_gt.compute_y_3_object(); + auto z_coord = m_gt.compute_z_3_object(); + auto vertex = m_gt.construct_vertex_3_object(); + + const Point_3& min_p = vertex(m_bbox, 0); + std::size_t i = (x_coord(p) - x_coord(min_p)) / m_spacing[0]; + std::size_t j = (y_coord(p) - y_coord(min_p)) / m_spacing[1]; + std::size_t k = (z_coord(p) - z_coord(min_p)) / m_spacing[2]; + + // @todo check this + if(i == xdim() - 1) + --i; + if(j == ydim() - 1) + --j; + if(k == zdim() - 1) + --k; + + return {i, j, k}; + } + + // Geometry +public: + /** + * \brief returns the geometric position of the grid vertex described by a set of indices. + * + * Depending on the value of the template parameter `cache_points`, positions might not be stored + * but calculated using the lowest corner of the bounding box and grid spacing. * * \param i the index in the `x` direction * \param j the index in the `y` direction * \param k the index in the `z` direction * - * \return the stored value - * * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` */ Point_3 point(const std::size_t i, const std::size_t j, const std::size_t k) const { - auto x_coord = m_gt.compute_x_3_object(); - auto y_coord = m_gt.compute_y_3_object(); - auto z_coord = m_gt.compute_z_3_object(); - auto point = m_gt.construct_point_3_object(); - auto vertex = m_gt.construct_vertex_3_object(); + typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object(); + typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object(); + typename Geom_traits::Construct_vertex_3 vertex = m_gt.construct_vertex_3_object(); const Point_3& min_p = vertex(m_bbox, 0); return point(x_coord(min_p) + i * m_spacing[0], y_coord(min_p) + j * m_spacing[1], z_coord(min_p) + k * m_spacing[2]); } - - /** - * \brief gets the scalar value stored at the grid vertex described by a set of indices. - * - * \param i the index in the `x` direction - * \param j the index in the `y` direction - * \param k the index in the `z` direction - * - * \return the stored value - * - * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` - */ - FT value(const std::size_t i, - const std::size_t j, - const std::size_t k) const - { - return m_values[linear_index(i, j, k)]; - } - - /** - * \brief gets the scalar value stored at the grid vertex described by a set of indices. - * - * \note This function can be used to set the value at a grid vertex. - * - * \param i the index in the `x` direction - * \param j the index in the `y` direction - * \param k the index in the `z` direction - * - * \return a reference to the stored value - * - * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` - */ - FT& value(const std::size_t i, - const std::size_t j, - const std::size_t k) - { - return m_values[linear_index(i, j, k)]; - } - -public: - /** - * \brief gets the gradient stored at the grid vertex described by a set of indices. - * - * \param i the index in the `x` direction - * \param j the index in the `y` direction - * \param k the index in the `z` direction - * - * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` - */ - const Vector_3& gradient(const std::size_t i, - const std::size_t j, - const std::size_t k) const - { - return m_gradients[linear_index(i, j, k)]; - } - - /** - * \brief gets the gradient stored at the grid vertex described by a set of indices. - * - * \note This function can be used to set the gradient at a grid vertex. - * - * \param i the index in the `x` direction - * \param j the index in the `y` direction - * \param k the index in the `z` direction - * - * \return a reference to the stored gradient - * - * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` - */ - Vector_3& gradient(const std::size_t i, - const std::size_t j, - const std::size_t k) - { - return m_gradients[linear_index(i, j, k)]; - } - -private: - std::size_t linear_index(const std::size_t i, - const std::size_t j, - const std::size_t k) const - { - CGAL_precondition(i < xdim() && j < ydim() && k < zdim()); - - // convert (i, j, k) into a linear index to access the scalar values / gradient vectors - return (k * ydim() + j) * xdim() + i; - } }; -template -Cartesian_grid_3:: -operator Image_3() const -{ - auto x_coord = m_gt.compute_x_3_object(); - auto y_coord = m_gt.compute_y_3_object(); - auto z_coord = m_gt.compute_z_3_object(); - auto vertex = m_gt.construct_vertex_3_object(); - - // select number type - WORD_KIND wordkind; - if(std::is_floating_point::value) // @fixme seems meaningless given that vx vy vz are doubles - 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; - - // get spacing - const double vx = CGAL::to_double(m_spacing[0]); - const double vy = CGAL::to_double(m_spacing[1]); - const double vz = CGAL::to_double(m_spacing[2]); - - // create image - _image* im = _createImage(xdim(), ydim(), zdim(), - 1, // vectorial dimension - vx, vy, vz, // voxel size - sizeof(FT), // image word size in bytes - 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 - const Point_3& min_p = vertex(m_bbox, 0); - im->tx = float(CGAL::to_double(x_coord(min_p))); - im->ty = float(CGAL::to_double(y_coord(min_p))); - im->tz = float(CGAL::to_double(z_coord(min_p))); - - // copy data - FT* data = static_cast(im->data); // @fixme what compatibility with non trivial FTs? - for(std::size_t x=0; x -bool write_OBJ(const std::string& filename, - const Cartesian_grid_3& grid, - const NamedParameters& np = parameters::default_values()) -{ - using Point_3 = typename GeomTraits::Point_3; - - auto x_coord = grid.geom_traits().compute_x_3_object(); - auto y_coord = grid.geom_traits().compute_y_3_object(); - auto z_coord = grid.geom_traits().compute_z_3_object(); - auto vertex = grid.geom_traits().construct_vertex_3_object(); - - std::ofstream out(filename); - set_ascii_mode(out); // obj is ASCII only - - set_stream_precision_from_NP(out, np); - - if(out.fail()) - return false; - - - // write vertices - for(std::size_t x=0; x + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domains_grp + * + * \cgalModels{IsosurfacingDomainWithGradient_3} + * + * \brief A domain that can be used as input in the %Dual Contouring algorithm. + * + * \details This class is essentially wrapper around the different bricks provided by its + * template parameters: `Partition` provides the spacial partitioning, `ValueField` and `GradientField` + * the values and gradients that define the isosurface. The optional template parameter + * `EdgeIntersectionOracle` gives control over the method used to computate edge-isosurface intersection points. + * + * \tparam Partition must be a model of `Partition_3` + * \tparam ValueField must be a model of `ValueField_3` + * \tparam GradientField must be a model of `GradientField_3` + * \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3` + * + * \sa `CGAL::Isosurfacing::dual_contouring()` + * \sa `CGAL::Isosurfacing::Marching_cubes_domain_3()` + */ +template +class Dual_contouring_domain_3 +#ifndef DOXYGEN_RUNNING + : public internal::Isosurfacing_domain_3 +#endif +{ +private: + using Base = internal::Isosurfacing_domain_3; + +public: + /** + * \brief constructs a domain that can be used with the %Dual Contouring algorithm. + * + * \param partition the space partitioning data structure + * \param values a continuous field of scalar values, defined over the bounding box of `partition` + * \param gradients a continuous field of normalized vectors, defined over the bounding box of `partition` + * \param intersection_oracle the oracle for edge-isosurface intersection computation + * + * \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects. + * As such, users must ensure that the lifetime of these objects exceeds that of the domain object. + */ + Dual_contouring_domain_3(const Partition& partition, + const ValueField& values, + const GradientField& gradients, + const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle()) + : Base(partition, values, gradients, intersection_oracle) + { } +}; + +/** + * \ingroup IS_Domains_grp + * + * \brief creates a new instance of a domain that can be used with the %Dual Contouring algorithm. + * + * \tparam Partition must be a model of `Partition_3` + * \tparam ValueField must be a model of `ValueField_3` + * \tparam GradientField must be a model of `GradientField_3` + * \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3` + * + * \param partition the space partitioning data structure + * \param values a continuous field of scalar values, defined over the bounding box of `partition` + * \param gradients a continuous field of normalized vectors, defined over the bounding box of `partition` + * \param intersection_oracle the oracle for edge-isosurface intersection computation + * + * \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects. + * As such, users must ensure that the lifetime of these objects exceeds that of the domain object. + */ +template +Dual_contouring_domain_3 +create_dual_contouring_domain_3(const Partition& partition, + const ValueField& values, + const GradientField& gradients, + const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle()) +{ + return { partition, values, gradients, intersection_oracle }; +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_DUAL_CONTOURING_DOMAIN_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h deleted file mode 100644 index 3c249ef9530..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_domain_3.h +++ /dev/null @@ -1,107 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H -#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H - -#include - -#include -#include -#include -#include -#include - -namespace CGAL { -namespace Isosurfacing { - -/** - * \ingroup IS_Domains_grp - * - * \cgalModels{IsosurfacingDomain_3,IsosurfacingDomainWithGradient_3} - * - * \brief A domain that represents an explicitly stored %Cartesian grid. - * - * \warning The domain keeps a pointer to the `grid` object, hence users must ensure that - * the lifetime of the `grid` object exceeds that of this object. - * - * \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter - * is a model of `IsosurfacingTraits_3`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` - * and implement `%Grid::GeomTraits::Vector_3 operator()(const %Grid::GeomTraits::Point_3& point) const`. - * - * \sa `CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain()` - */ -template - , typename Function = internal::Explicit_Cartesian_grid_function -#endif - > -class Explicit_Cartesian_grid_domain_3 -#ifndef DOXYGEN_RUNNING - : public internal::Isosurfacing_domain_3 -#endif -{ -private: - using Base = internal::Isosurfacing_domain_3; - -public: - /** - * \brief creates a domain that can be used as input for isosurfacing algorithms. - * - * \param grid the %Cartesian grid containing input data - * \param gradient a function giving the value of the gradient at each discretization point - */ - Explicit_Cartesian_grid_domain_3(const Grid& grid, - const Gradient& gradient = Gradient()) - : Base(Topology { grid.xdim(), grid.ydim(), grid.zdim() }, - Geometry { grid }, - Function { grid }, - gradient, - grid.geom_traits()) - { - } -}; - -/** - * \ingroup IS_Domains_grp - * - * \brief creates a domain that can be used as input for isosurfacing algorithms. - * - * \warning The domain keeps a pointer to the `grid` object, hence users must ensure that - * the lifetime of the `grid` object exceeds that of the object returned by this function. - * - * \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter - * is a model of `IsosurfacingTraits_3`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` - * and implement `%Grid::Geom_traits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - * - * \param grid the %Cartesian grid containing input data - * \param gradient a function giving the value of the gradient of the implicit function at each discretization point - */ -template -Explicit_Cartesian_grid_domain_3 -create_explicit_Cartesian_grid_domain(const Grid& grid, - const Gradient& gradient = Gradient()) -{ - return { grid, gradient }; -} - -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_DOMAIN_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h deleted file mode 100644 index 6f7b82b6830..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_Cartesian_grid_gradient_3.h +++ /dev/null @@ -1,138 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H -#define CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H - -#include - -#include - -namespace CGAL { -namespace Isosurfacing { - -/** - * \ingroup IS_Domains_grp - * - * \brief Class template for a gradient that is stored in a %Cartesian grid. - * - * \details The gradient at a query point is calculated by trilinear interpolation. - * - * \warning This class keeps a pointer to the `grid` object, hence users must ensure that - * the lifetime of the `grid` object exceeds that of this gradient object. - * - * \tparam Grid must be a `CGAL::Isosurfacing::Cartesian_grid_3` whose `GeomTraits` template parameter - * is a model of `IsosurfacingTraits_3`. -*/ -template // allow more than just Cartesian_grid_3 -class Explicit_Cartesian_grid_gradient_3 -{ -public: - using Geom_traits = typename Grid::Geom_traits; - using FT = typename Geom_traits::FT; - using Point_3 = typename Geom_traits::Point_3; - using Vector_3 = typename Geom_traits::Vector_3; - using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3; - -private: - const Grid& m_grid; - -public: - /** - * \brief creates a new instance of this gradient. - * - * \param grid the %Cartesian grid that stores the gradient - */ - Explicit_Cartesian_grid_gradient_3(const Grid& grid) - : m_grid(grid) - { } - - /** - * \brief evaluates the gradient at a point in space. - * - * \param p the position at which the gradient is computed - */ - Vector_3 operator()(const Point_3& p) const - { - typename Geom_traits::Compute_x_3 x_coord = m_grid.geom_traits().compute_x_3_object(); - typename Geom_traits::Compute_y_3 y_coord = m_grid.geom_traits().compute_y_3_object(); - typename Geom_traits::Compute_z_3 z_coord = m_grid.geom_traits().compute_z_3_object(); - typename Geom_traits::Construct_vector_3 vector = m_grid.geom_traits().construct_vector_3_object(); - typename Geom_traits::Construct_vertex_3 vertex = m_grid.geom_traits().construct_vertex_3_object(); - - // trilinear interpolation of stored gradients - const Iso_cuboid_3& bbox = m_grid.bbox(); - const std::array& spacing = m_grid.spacing(); - - // calculate min index including border case - const Point_3& min_p = vertex(bbox, 0); - std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0]; - std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1]; - std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2]; - - // @todo check this - if(i == m_grid.xdim() - 1) - --i; - if(j == m_grid.ydim() - 1) - --j; - if(k == m_grid.zdim() - 1) - --k; - - // calculate coordinates of min index - const FT min_x = i * spacing[0] + x_coord(min_p); - const FT min_y = j * spacing[1] + y_coord(min_p); - const FT min_z = k * spacing[2] + z_coord(min_p); - - // interpolation factors between 0 and 1 - const FT f_i = (x_coord(p) - min_x) / spacing[0]; - const FT f_j = (y_coord(p) - min_y) / spacing[1]; - const FT f_k = (z_coord(p) - min_z) / spacing[2]; - - // read the gradient at all 8 corner points - const Vector_3& g000 = m_grid.gradient(i + 0, j + 0, k + 0); - const Vector_3& g001 = m_grid.gradient(i + 0, j + 0, k + 1); - const Vector_3& g010 = m_grid.gradient(i + 0, j + 1, k + 0); - const Vector_3& g011 = m_grid.gradient(i + 0, j + 1, k + 1); - const Vector_3& g100 = m_grid.gradient(i + 1, j + 0, k + 0); - const Vector_3& g101 = m_grid.gradient(i + 1, j + 0, k + 1); - const Vector_3& g110 = m_grid.gradient(i + 1, j + 1, k + 0); - const Vector_3& g111 = m_grid.gradient(i + 1, j + 1, k + 1); - - // interpolate along all axes by weighting the corner points - const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k); - const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k; - const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k); - const FT lambda011 = (FT(1) - f_i) * f_j * f_k; - const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k); - const FT lambda101 = f_i * (FT(1) - f_j) * f_k; - const FT lambda110 = f_i * f_j * (FT(1) - f_k); - const FT lambda111 = f_i * f_j * f_k; - - // add weighted corners - return vector( x_coord(g000) * lambda000 + x_coord(g001) * lambda001 + - x_coord(g010) * lambda010 + x_coord(g011) * lambda011 + - x_coord(g100) * lambda100 + x_coord(g101) * lambda101 + - x_coord(g110) * lambda110 + x_coord(g111) * lambda111, - y_coord(g000) * lambda000 + y_coord(g001) * lambda001 + - y_coord(g010) * lambda010 + y_coord(g011) * lambda011 + - y_coord(g100) * lambda100 + y_coord(g101) * lambda101 + - y_coord(g110) * lambda110 + y_coord(g111) * lambda111, - z_coord(g000) * lambda000 + z_coord(g001) * lambda001 + - z_coord(g010) * lambda010 + z_coord(g011) * lambda011 + - z_coord(g100) * lambda100 + z_coord(g101) * lambda101 + - z_coord(g110) * lambda110 + z_coord(g111) * lambda111 ); - } -}; - -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_EXPLICIT_CARTESIAN_GRID_GRADIENT_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Finite_difference_gradient_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Finite_difference_gradient_3.h index de14ac26ae2..7dfdedd2974 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Finite_difference_gradient_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Finite_difference_gradient_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,24 +8,29 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Julian Stahl +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H #define CGAL_ISOSURFACING_3_FINITE_DIFFERENCE_GRADIENT_3_H #include +#include + namespace CGAL { namespace Isosurfacing { /** - * \ingroup IS_Domain_helpers_grp + * \ingroup IS_Fields_grp + * + * \cgalModels{GradientField_3} * * \brief Class template for a gradient that is calculated using finite differences. * - * \details This gradient function evaluates an implicit value function at six points - * that are a given distance `delta` away from the queried point along the %Cartesian axes. + * \details This gradient function evaluates a value function at six points that are + * a given distance `delta` away from the queried point along the %Cartesian axes. * - * \tparam GeomTraits must be a model of `Kernel`. + * \tparam GeomTraits must be a model of `IsosurfacingTraits_3`. */ template class Finite_difference_gradient_3 @@ -44,12 +49,11 @@ private: public: /** - * \brief creates a new instance of this gradient. + * \brief creates a new instance of this gradient class. * - * \tparam ValueFunction the type of the implicit function. It must be a model of `CopyConstructible` - * and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. + * \tparam ValueFunction must be a model of `ValueFunction_3`. * - * \param func the implicit function giving the value of the implicit function at each discretization point + * \param function the function giving the scalar value at each point * \param delta the distance for calculating the finite differences * \param gt the geometric traits class */ @@ -64,7 +68,7 @@ public: { } /** - * \brief evaluates the gradient at a point in space. + * \brief evaluates the gradient at a point in 3D space. * * \param p the position at which the gradient is computed */ diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Gradient_function_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Gradient_function_3.h new file mode 100644 index 00000000000..6f8cc5d290a --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Gradient_function_3.h @@ -0,0 +1,90 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H +#define CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H + +#include + +#include + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Fields_grp + * + * \cgalModels{GradientField_3} + * + * \brief The class `Gradient_function_3` represents a field of vectors computed + * using a user-provided unary function. + * + * \tparam Partition must be a model of `Partition_3` + * + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3` + */ +template +class Gradient_function_3 +{ +public: + using Geom_traits = typename Partition::Geom_traits; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; + + using PT = partition_traits; + using Vertex_descriptor = typename PT::Vertex_descriptor; + +private: + std::function m_fn; + const Partition& m_partition; + +public: + /** + * \brief constructs a field of gradients using a gradient function and a partition. + * + * \tparam Function must provide the following function signature: + * `Vector_3 operator()(const %Point_3&) const` + * + * \param fn the function providing gradients + * \param partition the space partitioning data structure + */ + template + Gradient_function_3(const Function& fn, + const Partition& partition) + : m_fn{fn}, + m_partition{partition} + { } + +public: + /** + * \brief evaluates the function at the point `p`. + */ + Vector_3 operator()(const Point_3& p) const + { + return m_fn(p); + } + + /** + * \brief evaluates the function at the vertex `v`. + */ + const Vector_3& operator()(const Vertex_descriptor& v) const + { + return this->operator()(PT::point(v, m_partition)); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_GRADIENT_FUNCTION_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h new file mode 100644 index 00000000000..f1bc964004d --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/IO/Image_3.h @@ -0,0 +1,170 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_IO_IMAGE_3_H +#define CGAL_ISOSURFACING_3_IO_IMAGE_3_H + +#include + +#include +#include + +#include + +namespace CGAL { +namespace Isosurfacing { +namespace IO { + +/** + * \ingroup IS_IO_functions_grp + * + * \brief creates a grid from a `CGAL::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. + * + * \tparam K must be a model of `IsosurfacingTraits_3` + * + * \param image the image providing the data + * \param k the traits + */ +template +std::pair, + Interpolated_discrete_values_3 > > +read_Image_3(const CGAL::Image_3& image, + const K& k = K()) +{ + using Grid = Cartesian_grid_3; + using Values = Interpolated_discrete_values_3; + + using FT = typename K::FT; + using Point_3 = typename K::Point_3; + using Iso_cuboid_3 = typename K::Iso_cuboid_3; + + auto point = k.construct_point_3_object(); + auto iso_cuboid = k.construct_iso_cuboid_3_object(); + + Iso_cuboid_3 bbox; + + // 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 = iso_cuboid(point(image.tx(), image.ty(), image.tz()), + point(max_x, max_y, max_z)); + + // get spacing + // std::array spacing = make_array(image.vx(), image.vy(), image.vz()); + + // get sizes + std::array sizes; + sizes[0] = image.xdim(); + sizes[1] = image.ydim(); + sizes[2] = image.zdim(); + + Grid grid { bbox, sizes[0], sizes[1], sizes[2], k }; + + Values values { grid }; + + // copy values + for(std::size_t x=0; x` with `GeomTraits` + * a model of `IsosurfacingTraits_3` + * + * \param grid the space partitioning data structure + * \param values the field of values + */ +template +CGAL::Image_3 write_Image_3(const Grid& grid, + const Values& values) +{ + using Geom_traits = typename Grid::Geom_traits; + + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + + const Geom_traits& gt = grid.geom_traits(); + typename Geom_traits::Compute_x_3 x_coord = gt.compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = gt.compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = gt.compute_z_3_object(); + typename Geom_traits::Construct_vertex_3 vertex = gt.construct_vertex_3_object(); + + // select number type + WORD_KIND wordkind; + if(std::is_floating_point::value) // @fixme seems meaningless given that vx vy vz are doubles + 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; + + // get spacing + const double vx = CGAL::to_double(grid.spacing()[0]); + const double vy = CGAL::to_double(grid.spacing()[1]); + const double vz = CGAL::to_double(grid.spacing()[2]); + + // create image + _image* im = _createImage(grid.xdim(), grid.ydim(), grid.zdim(), + 1, // vectorial dimension + vx, vy, vz, // voxel size + sizeof(FT), // image word size in bytes + 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 + const Point_3& min_p = vertex(grid.bbox(), 0); + im->tx = float(CGAL::to_double(x_coord(min_p))); + im->ty = float(CGAL::to_double(y_coord(min_p))); + im->tz = float(CGAL::to_double(z_coord(min_p))); + + // copy data + FT* data = static_cast(im->data); // @fixme what compatibility with non trivial FTs? + for(std::size_t x=0; x + +#include + +#include +#include + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +template +class Cartesian_grid_3; + +namespace IO { + +template +bool write_OBJ(std::ostream& out, + const Cartesian_grid_3& grid, + const NamedParameters& np = parameters::default_values()) +{ + using Point_3 = typename GeomTraits::Point_3; + + auto x_coord = grid.geom_traits().compute_x_3_object(); + auto y_coord = grid.geom_traits().compute_y_3_object(); + auto z_coord = grid.geom_traits().compute_z_3_object(); + auto vertex = grid.geom_traits().construct_vertex_3_object(); + + set_ascii_mode(out); // obj is ASCII only + + set_stream_precision_from_NP(out, np); + + if(out.fail()) + return false; + + // write vertices + for(std::size_t x=0; x +bool write_OBJ(const std::string& fname, + const Cartesian_grid_3& grid, + const NamedParameters& np = parameters::default_values()) +{ + std::ofstream os(fname); + CGAL::IO::set_mode(os, CGAL::IO::ASCII); + + return write_OBJ(os, grid, np); +} + +} // namespace IO +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_IMAGE_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h deleted file mode 100644 index 95d2f92f34d..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_Cartesian_grid_domain_3.h +++ /dev/null @@ -1,172 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H -#define CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H - -#include - -#include -#include -#include -#include -#include - -#include - -namespace CGAL { -namespace Isosurfacing { - -/** - * \ingroup IS_Domains_grp - * - * \cgalModels{IsosurfacingDomain_3,IsosurfacingDomainWithGradient_3} - * - * \brief A domain that represents a %Cartesian grid that discretizes an implicit function. - * - * \tparam GeomTraits must be a model of `IsosurfacingTraits_3`. - * \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` - * and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement - * `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - * - * \sa `CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain()` - */ -template - , typename Function = internal::Implicit_function_with_geometry -#endif - > -class Implicit_Cartesian_grid_domain_3 -#ifndef DOXYGEN_RUNNING - : public internal::Isosurfacing_domain_3 -#endif -{ -private: - using Base = internal::Isosurfacing_domain_3; - - using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; - using Vector_3 = typename GeomTraits::Vector_3; - - Base construct_domain(const typename GeomTraits::Iso_cuboid_3& bbox, - const typename GeomTraits::Vector_3& spacing, - const ImplicitFunction& point_function, - const Gradient& gradient, - const GeomTraits& gt) - { - auto x_coord = gt.compute_x_3_object(); - auto y_coord = gt.compute_y_3_object(); - auto z_coord = gt.compute_z_3_object(); - auto vertex = gt.construct_vertex_3_object(); - auto vector = gt.construct_vector_3_object(); - - const Point_3& min_p = vertex(bbox, 0); - const Point_3& max_p = vertex(bbox, 7); - const FT x_span = x_coord(max_p) - x_coord(min_p); - const FT y_span = y_coord(max_p) - y_coord(min_p); - const FT z_span = z_coord(max_p) - z_coord(min_p); - - const std::size_t x_dim = std::ceil(x_span / x_coord(spacing)) + 1; - const std::size_t y_dim = std::ceil(y_span / y_coord(spacing)) + 1; - const std::size_t z_dim = std::ceil(z_span / z_coord(spacing)) + 1; - - Topology topology { x_dim, y_dim, z_dim }; - - const Vector_3 offset = vector(x_coord(min_p), y_coord(min_p), z_coord(min_p)); - Geometry geometry { offset, spacing }; - - Function func { geometry, point_function }; - - return { topology, geometry, func, gradient, gt }; - } - -public: - /** - * \brief creates a domain from an implicit function. - * - * \details The implicit function is evaluated at the vertices of the virtual grid - * defined by the bounding box and the spacing value. By not storing any function values explicitely, - * less overall memory is required in comparison to an `Explicit_Cartesian_grid_domain_3`. - * - * \tparam GeomTraits must be a model of `IsosurfacingTraits_3`. - * \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` - * and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement - * `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - * - * \param bbox an axis-aligned box that specifies the dimensions of the implicit function's domain - * \param spacing the distance between discretization points - * \param point_function the implicit function giving the value of the implicit function at each discretization point - * \param gradient a function giving the value of the gradient of the implicit function at each discretization point - * \param gt an instance of geometric traits - * - * \pre `spacing != CGAL::NULL_VECTOR` - */ - Implicit_Cartesian_grid_domain_3(const typename GeomTraits::Iso_cuboid_3& bbox, - const typename GeomTraits::Vector_3& spacing, - const ImplicitFunction& point_function, - const Gradient& gradient = Gradient(), - const GeomTraits& gt = GeomTraits()) - : Base(construct_domain(bbox, spacing, point_function, gradient, gt)) - { - } -}; - -/** - * \ingroup IS_Domains_grp - * - * \brief creates a domain from an implicit function that can be used as input for isosurfacing algorithms. - * - * \details The implicit function is evaluated at the vertices of the virtual grid - * defined by the bounding box and the spacing value. By not storing any function values explicitely, - * less overall memory is required in comparison to an `Explicit_Cartesian_grid_domain_3`. - * - * \tparam GeomTraits must be a model of `IsosurfacingTraits_3`. - * \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` - * and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement - * `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - * - * \param bbox an axis-aligned box that specifies the dimensions of the implicit function's domain - * \param spacing the distance between discretization points - * \param point_function the implicit function giving the value of the implicit function at each discretization point - * \param gradient a function giving the value of the gradient of the implicit function at each discretization point - * \param gt an instance of geometric traits - * - * \return a new instance of `CGAL::Isosurfacing::Implicit_Cartesian_grid_domain_3` - * - * \pre `spacing != CGAL::NULL_VECTOR` - */ -template -Implicit_Cartesian_grid_domain_3 -create_implicit_Cartesian_grid_domain(const typename GeomTraits::Iso_cuboid_3& bbox, - const typename GeomTraits::Vector_3& spacing, - const ImplicitFunction& point_function, - const Gradient& gradient = Gradient(), - const GeomTraits& gt = GeomTraits()) -{ - return { bbox, spacing, point_function, gradient, gt }; -} - -// @todo add an undocumented convenience overload with Vector_3 to match CGAL kernels -// without having to provide the kernel in the call like f(...) - -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_IMPLICIT_CARTESIAN_GRID_DOMAIN_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h deleted file mode 100644 index e7e9a0273ec..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h +++ /dev/null @@ -1,108 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H -#define CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H - -#include - -#include -#include -#include -#include -#include -#include - -namespace CGAL { -namespace Isosurfacing { - -/* - * \ingroup IS_Domains_grp - * - * \cgalModels{IsosurfacingDomainWithGradient_3} - * - * \brief A domain that represents an octree that discretizes an implicit function. - * - * \warning The domain keeps a pointer to the `octree` object, hence users must ensure that - * the lifetime of the `octree` object exceeds that of this object. - * - * \tparam Octree must be a `CGAL::Octree`. - * \tparam ImplicitFunction the type of the implicit function. It must be a model of CopyConstructible implement - * `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement - * `GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - */ -template - , typename Geometry = internal::Octree_geometry - , typename Function = internal::Implicit_function_with_geometry -#endif - > -class Implicit_octree_domain -#ifndef DOXYGEN_RUNNING - : public internal::Isosurfacing_domain_3 -#endif -{ - using Base = internal::Isosurfacing_domain_3; - -public: - Implicit_octree_domain(const Octree& octree, - const ImplicitFunction& point_function, - const Gradient& gradient = Gradient()) - : Base(Topology { octree }, - Geometry { octree }, - Function { Geometry { octree }, point_function }, - gradient, - octree.geom_traits()) - { - } -}; - -/* - * \ingroup IS_Domains_grp - * - * \brief creates a domain from an octree that can be used as input for isosurfacing algorithms. - * - * \warning The domain keeps a pointer to the `octree` object, hence users must ensure that - * the lifetime of the `octree` object exceeds that of the object returned by this function. - * - * \tparam Octree must be a `CGAL::Octree`. - * \tparam ImplicitFunction the type of the implicit function. It must be a model of `CopyConstructible` - * and implement `GeomTraits::FT operator()(const GeomTraits::Point_3& point) const`. - * \tparam Gradient the type of the gradient functor. It must be a model of `CopyConstructible` and implement - * `Octree::GeomTraits::Vector_3 operator()(const GeomTraits::Point_3& point) const`. - * - * \param octree an octree - * \param point_function the function with a point as argument - * \param gradient a function that describes the gradient of the data - * - * \return a new `CGAL::Implicit_octree_domain` - */ -template -Implicit_octree_domain -create_implicit_octree_domain(const Octree& octree, - const ImplicitFunction& point_function, - const Gradient& gradient = Gradient()) -{ - return { octree, point_function, gradient }; -} - -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_IMPLICIT_OCTREE_DOMAIN_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h new file mode 100644 index 00000000000..a3fd48e02d6 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h @@ -0,0 +1,124 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H +#define CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H + +#include + +#include + +#include +#include + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{ValueField_3} + * + * \brief Class template for a gradient that is calculated using discrete values and trilinear interpolation. + * + * \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3`, with `GeomTraits` a model of `IsosurfacingTraits_3` + * \tparam InterpolationScheme must be a model of `InterpolationScheme_3` + */ +template > +class Interpolated_discrete_gradients_3 +{ + using Geom_traits = typename Grid::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; + + using Vertex_descriptor = typename partition_traits::Vertex_descriptor; + +private: + const Grid& m_grid; + const InterpolationScheme m_interpolation; + + std::vector m_gradients; + +public: + Interpolated_discrete_gradients_3(const Grid& grid, + const InterpolationScheme& interpolation = InterpolationScheme()) + : m_grid{grid}, + m_interpolation{interpolation} + { + // pre-allocate memory + const std::size_t nv = grid.xdim() * grid.ydim() * grid.zdim(); + m_gradients.resize(nv); + } + + // computes and stores gradients at the vertices of the grid + // \tparam must be ValueField a model of `ValueField_3` + // \param values a field of values whose gradient are being computed + template + void compute_discrete_gradients(const ValueField& values) + { + // @todo + } + +public: + /** + * \brief gets the gradient stored at the grid vertex described by a set of indices. + * + * \note This function can be used to set the gradient at a grid vertex. + * + * \param i the index in the `x` direction + * \param j the index in the `y` direction + * \param k the index in the `z` direction + * + * \return a reference to the stored gradient + * + * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` + */ + Vector_3& operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) + { + return m_gradients[m_grid.linear_index(i, j, k)]; + } + + /** + * \brief gets the gradient stored at the grid vertex described by a set of indices. + * + * \param i the index in the `x` direction + * \param j the index in the `y` direction + * \param k the index in the `z` direction + * + * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` + */ + const Vector_3& operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) const + { + return m_gradients[m_grid.linear_index(i, j, k)]; + } + + /*! + * returns the gradient at a given point `p`. + */ + Vector_3 operator()(const Point_3& p) const + { + return m_interpolation.interpolate_gradients(p, m_grid, m_gradients); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_GRADIENTS_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h new file mode 100644 index 00000000000..ae9b117aca1 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h @@ -0,0 +1,124 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H +#define CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H + +#include + +#include + +#include + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{ValueField_3} + * + * \brief Class template for a field of values that are calculated using discrete values and trilinear interpolation. + * + * \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3`, with `GeomTraits` a model of `IsosurfacingTraits_3` + * \tparam InterpolationScheme must be a model of `InterpolationScheme_3` + */ +template > +class Interpolated_discrete_values_3 +{ + using Geom_traits = typename Grid::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + + using Vertex_descriptor = typename partition_traits::Vertex_descriptor; + +private: + const Grid& m_grid; + const InterpolationScheme m_interpolation; + + std::vector m_values; + +public: + Interpolated_discrete_values_3(const Grid& grid, + const InterpolationScheme& interpolation = InterpolationScheme()) + : m_grid{grid}, + m_interpolation{interpolation} + { + // pre-allocate memory + const std::size_t nv = grid.xdim() * grid.ydim() * grid.zdim(); + m_values.resize(nv); + } + +public: + /** + * \brief gets the scalar value stored at the grid vertex described by a set of indices. + * + * \note This function can be used to set the value at a grid vertex. + * + * \param i the index in the `x` direction + * \param j the index in the `y` direction + * \param k the index in the `z` direction + * + * \return a reference to the stored value + * + * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` + */ + FT& operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) + { + return m_values[m_grid.linear_index(i, j, k)]; + } + + /** + * \brief gets the scalar value stored at the grid vertex described by a set of indices. + * + * \param i the index in the `x` direction + * \param j the index in the `y` direction + * \param k the index in the `z` direction + * + * \return the stored value + * + * \pre `i < xdim()` and `j < ydim()` and `k < zdim()` + */ + FT operator()(const std::size_t i, + const std::size_t j, + const std::size_t k) const + { + return m_values[m_grid.linear_index(i, j, k)]; + } + + /*! + * returns the value at vertex `v`. + */ + FT operator()(const Vertex_descriptor& v) const + { + return this->operator()(v[0], v[1], v[2]); + } + + /*! + * returns the value at point `p`. + */ + FT operator()(const Point_3& p) const + { + return m_interpolation.interpolate_values(p, m_grid, m_values); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERPOLATED_DISCRETE_VALUES_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Marching_cubes_domain_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Marching_cubes_domain_3.h new file mode 100644 index 00000000000..f8672cd200b --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Marching_cubes_domain_3.h @@ -0,0 +1,102 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H +#define CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H + +#include + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domains_grp + * + * \cgalModels{IsosurfacingDomain_3} + * + * \brief A domain that can be used with the Marching Cubes algorithm. + * + * \details This class is essentially wrapper around the different bricks provided by its + * template parameters: `Partition` provides the spacial partitioning, `ValueField` + * the values that define the isosurface. The optional template parameter + * `EdgeIntersectionOracle` gives control over the method used to computate edge-isosurface intersection points. + * + * \tparam Partition must be a model of `Partition_3` + * \tparam ValueField must be a model of `ValueField_3` + * \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3` + * + * \sa `CGAL::Isosurfacing::marching_cubes_3()` + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3` + */ +template +class Marching_cubes_domain_3 +#ifndef DOXYGEN_RUNNING + : public internal::Isosurfacing_domain_3 +#endif +{ +private: + using Base = internal::Isosurfacing_domain_3; + +public: + /** + * \brief constructs a domain that can be used with the Marching Cubes algorithm. + * + * \param partition the space partitioning data structure + * \param values a continuous field of scalar values, defined over the bounding box of `partition` + * \param intersection_oracle the oracle for edge-isosurface intersection computation + * + * \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects. + * As such, users must ensure that the lifetime of these objects exceeds that of the domain object. + */ + Marching_cubes_domain_3(const Partition& partition, + const ValueField& values, + const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle()) + : Base(partition, values, intersection_oracle) + { } +}; + +/** + * \ingroup IS_Domains_grp + * + * \brief creates a new instance of a domain that can be used with the Marching Cubes algorithm. + * + * \tparam Partition must be a model of `Partition_3` + * \tparam ValueField must be a model of `ValueField_3` + * \tparam EdgeIntersectionOracle must be a model of `EdgeIntersectionOracle_3` + * + * \param partition the space partitioning data structure + * \param values a continuous field of scalar values, defined over the bounding box of `partition` + * \param intersection_oracle the oracle for edge-isosurface intersection computation + * + * \warning the domain class keeps a reference to the `partition`, `values` and `gradients` objects. + * As such, users must ensure that the lifetime of these objects exceeds that of the domain object. + */ +template +Marching_cubes_domain_3 +create_marching_cubes_domain_3(const Partition& partition, + const ValueField& values, + const EdgeIntersectionOracle& intersection_oracle = EdgeIntersectionOracle()) +{ + return { partition, values, intersection_oracle }; +} + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_MARCHING_CUBES_DOMAIN_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Value_function_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Value_function_3.h new file mode 100644 index 00000000000..e82b3dbafdb --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Value_function_3.h @@ -0,0 +1,91 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H +#define CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H + +#include + +#include + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{ValueField_3} + * + * \brief The class `Value_function_3` represents a field of scalars computed + * using a user-provided unary function. + * + * \tparam Partition must be a model of `Partition_3` + * + * \sa `CGAL::Isosurfacing::Marching_cubes_domain_3` + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3` + */ +template +class Value_function_3 +{ +public: + using Geom_traits = typename Partition::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + + using PT = partition_traits; + using Vertex_descriptor = typename PT::Vertex_descriptor; + +private: + std::function m_fn; + const Partition& m_partition; + +public: + /** + * \brief constructs a field of values using a value function and a partition. + * + * \tparam Function must provide the following function signature: + * `FT operator()(const %Point_3&) const` + * + * \param fn the function providing values + * \param partition the space partitioning data structure + */ + template + Value_function_3(const Function& fn, + const Partition& partition) + : m_fn{fn}, + m_partition{partition} + { } + +public: + /** + * \brief evaluates the function at the point `p` + */ + FT operator()(const Point_3& p) const + { + return m_fn(p); + } + + /** + * \brief evaluates the function at the vertex `v` + */ + FT operator()(const Vertex_descriptor& v) const + { + return this->operator()(PT::point(v, m_partition)); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_VALUE_FUNCTION_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Zero_gradient.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Zero_gradient.h deleted file mode 100644 index 645814db1b7..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Zero_gradient.h +++ /dev/null @@ -1,46 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_ZERO_GRADIENT_H -#define CGAL_ISOSURFACING_3_ZERO_GRADIENT_H - -#include - -#include - -namespace CGAL { -namespace Isosurfacing { - -/** - * \ingroup IS_Domain_helpers_grp - * - * \brief Class template for a gradient that equals zero everywhere. - * - * \tparam P the point type - * - * \details This gradient function can be useful when using Marching Cubes, which does not require the domain to have a gradient. - */ -struct Zero_gradient -{ - /** - * \return the null vector - */ - template - CGAL::Null_vector operator()(const P&) const - { - return CGAL::NULL_VECTOR; - } -}; - -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_ZERO_GRADIENT_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/dual_contouring_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/dual_contouring_3.h index 34df686ac9b..ad6d508ff1b 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/dual_contouring_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/dual_contouring_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -9,6 +9,7 @@ // // Author(s) : Julian Stahl // Daniel Zint +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_DUAL_CONTOURING_3_H #define CGAL_ISOSURFACING_3_DUAL_CONTOURING_3_H @@ -17,7 +18,6 @@ #include -#include #include namespace CGAL { @@ -26,7 +26,10 @@ namespace Isosurfacing { /** * \ingroup IS_Methods_grp * - * \brief creates an indexed face set that represents an isosurface generated by the %Dual Contouring algorithm. + * \brief creates a polygon soup that discretizes an isosurface using the %Dual Contouring algorithm. + * + * \details The point placement strategy within each cell of the space partition is based on + * Quadric Error Metrics ("QEM", or "QEF" in %Dual Contouring-related works). * * \tparam ConcurrencyTag enables sequential versus parallel algorithm. * Possible values are `Sequential_tag`, `Parallel_if_available_tag`, or `Parallel_tag`. @@ -36,63 +39,54 @@ namespace Isosurfacing { * \tparam PolygonRange must be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence` * whose value type is itself a model of the concepts `RandomAccessContainer` * and `BackInsertionSequence` whose value type is `std::size_t`. + * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * - * \param domain the domain providing input data and its topology - * \param isovalue value of the isosurface - * \param points points of the polygons in the created indexed face set + * \param domain the domain providing the spacial partition and the values and gradient data + * \param isovalue the value defining the isosurface + * \param points the points of the polygons in the created polygon soup * \param polygons each element in the vector describes a polygon using the indices of the points in `points` + * \param np optional \ref bgl_namedparameters "Named Parameters" described below + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{constrain_to_cell} + * \cgalParamDescription{whether to constrain the vertex position to the geometrical space of its cell} + * \cgalParamType{Boolean} + * \cgalParamDefault{`false`} + * \cgalParamExtra{Constraining the vertex to its dual cell guarantees that the resulting + * surface is a without self-intersections (non-manifoldness aside). Oppositely, + * an unconstrained positioning strategy might produce better looking surfaces + * near sharp features (ridges, corners), at the cost of possible self-intersections.} + * \cgalParamNEnd + * + * \cgalParamNBegin{do_not_triangulate_faces} + * \cgalParamDescription{If `true`, the output will contain quadrilaterals. + * If `false`, the output will contain triangles.} + * \cgalParamType{Boolean} + * \cgalParamDefault{`false` (faces are triangulated)} + * \cgalParamExtra{Triangulating faces is done by inserting the intersection between an edge and + * the isosurface, and linking it to the dual points of the cells incident to the edge. + * If `constrain_to_cell` is set to `false`, triangulation faces can result in additional + * self-intersections. An alternative that has worse approximation but is less likely + * to produce self-intersections is to use the function + * `CGAL::Polygon_mesh_processing::triangulate_faces()`.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` */ -#ifdef DOXYGEN_RUNNING // Do not document Positioning -template -void dual_contouring(const Domain& domain, - const typename Domain::Geom_traits::FT isovalue, - PointRange& points, - PolygonRange& polygons) -#else template > + typename NamedParameters = parameters::Default_named_parameters> void dual_contouring(const Domain& domain, const typename Domain::Geom_traits::FT isovalue, PointRange& points, PolygonRange& polygons, - const Positioning& positioning = Positioning()) -#endif + const NamedParameters& np = parameters::default_values()) { - // create vertices in each relevant cell - internal::Dual_contouring_vertex_positioning pos_func(domain, isovalue, positioning); - domain.template for_each_cell(pos_func); - - // connect vertices around an edge to form a face - internal::Dual_contouring_face_generation face_generation(domain, isovalue); - domain.template for_each_edge(face_generation); - - // copy vertices to point range - CGAL::internal::resize(points, 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 not stored into 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 degenerate faces - if(vertex_ids.size() > 2) - polygons.push_back(vertex_ids); - } + internal::Dual_contourer contourer; + contourer(domain, isovalue, points, polygons, np); } } // namespace Isosurfacing diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/edge_intersection_oracles_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/edge_intersection_oracles_3.h new file mode 100644 index 00000000000..bfd6ed04d5c --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/edge_intersection_oracles_3.h @@ -0,0 +1,212 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl +// Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H +#define CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +/** + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{EdgeIntersectionOracle_3} + * + * \brief The class `Dichotomy_edge_intersection` uses a dichotomy to find the intersection point + * between an edge and the isosurface. + * + * This class is for example suitable to be used as the `EdgeIntersectionOracle` template parameter of isosurfacing + * domain classes when values are computed using an implicit function. + * It is however not optimal when the values are interpolated from discrete values + * since the intersection can be computed analytically in this case. + * + * \sa `CGAL::Isosurfacing::Linear_interpolation_edge_intersection` + * \sa `CGAL::Isosurfacing::Marching_cubes_domain_3` + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3` + * \sa `CGAL::Isosurfacing::Value_field_3` + */ +struct Dichotomy_edge_intersection +{ + /*! + * \brief computes the intersection point between an edge and the isosurface. + * + * \tparam Domain must be a model of `IsosurfacingDomain_3` + * + * \param p_0 the geometric position of the first vertex of the edge + * \param p_1 the geometric position of the second vertex of the edge + * \param val_0 the value at the first vertex of the edge + * \param val_1 the value at the second vertex of the edge + * \param domain the isosurfacing domain + * \param isovalue the isovalue defining the isosurfacing with which we seek an intersection + * \param p the intersection point, if it exists + * + * \return `true` if the intersection point exists, `false` otherwise + */ + template // == Isosurfacing_domain_3 or similar + bool operator()(const typename Domain::Geom_traits::Point_3& p_0, + const typename Domain::Geom_traits::Point_3& p_1, + const typename Domain::Geom_traits::FT val_0, + const typename Domain::Geom_traits::FT val_1, + const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + typename Domain::Geom_traits::Point_3& p) const + { + using FT = typename Domain::Geom_traits::FT; + using Point_3 = typename Domain::Geom_traits::Point_3; + + using Vertex_descriptor = typename Domain::Vertex_descriptor; + using Edge_descriptor = typename Domain::Edge_descriptor; + + auto x_coord = domain.geom_traits().compute_x_3_object(); + auto y_coord = domain.geom_traits().compute_y_3_object(); + auto z_coord = domain.geom_traits().compute_z_3_object(); + auto point = domain.geom_traits().construct_point_3_object(); + + if((val_0 <= isovalue) == (val_1 <= isovalue)) + return false; + + Point_3 pl = p_0; + Point_3 pr = p_1; + FT val_l = val_0; + FT val_r = val_1; + + // @todo max iter choice + unsigned int dichotomy_iterations = 10, iter = 0; + + do + { + p = point((x_coord(pl) + x_coord(pr)) / FT(2), + (y_coord(pl) + y_coord(pr)) / FT(2), + (z_coord(pl) + z_coord(pr)) / FT(2)); + + const FT val_p = domain.value(p); + + if((val_l <= isovalue) != (val_p <= isovalue)) + { + pr = p; + val_r = val_p; + } + else + { + pl = p; + val_l = val_p; + } + + // @todo something like 1e-10 * isovalue? + // see also https://stats.stackexchange.com/questions/86708/how-to-calculate-relative-error-when-the-true-value-is-zero + if(is_zero(val_l - val_r)) + break; + } + while(++iter < dichotomy_iterations); + + + return true; + } +}; + +/** + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{EdgeIntersectionOracle_3} + * + * \brief The class `Linear_interpolation_edge_intersection` uses linear interpolation + * to find the intersection point between an edge and the isosurface. + * + * This class is for example suitable when interpolated discrete values are being used. + * + * \sa `CGAL::Isosurfacing::Dichotomy_edge_intersection` + * \sa `CGAL::Isosurfacing::Marching_cubes_domain_3` + * \sa `CGAL::Isosurfacing::Dual_contouring_domain_3` + * \sa `CGAL::Isosurfacing::Interpolated_discrete_values_3` + */ +struct Linear_interpolation_edge_intersection +{ + /*! + * \brief computes the intersection point between an edge and the isosurface. + * + * \tparam Domain must be a model of `IsosurfacingDomain_3` + * + * \param p_0 the geometric position of the first vertex of the edge + * \param p_1 the geometric position of the second vertex of the edge + * \param val_0 the value at the first vertex of the edge + * \param val_1 the value at the second vertex of the edge + * \param domain the isosurfacing domain + * \param isovalue the isovalue defining the isosurfacing with which we seek an intersection + * \param p the intersection point, if it exists + * + * \return `true` if the intersection point exists, `false` otherwise + */ + template + bool operator()(const typename Domain::Geom_traits::Point_3& p_0, + const typename Domain::Geom_traits::Point_3& p_1, + const typename Domain::Geom_traits::FT val_0, + const typename Domain::Geom_traits::FT val_1, + const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + typename Domain::Geom_traits::Point_3& p) const + { + using FT = typename Domain::Geom_traits::FT; + using Point_3 = typename Domain::Geom_traits::Point_3; + + using Vertex_descriptor = typename Domain::Vertex_descriptor; + using Edge_descriptor = typename Domain::Edge_descriptor; + + auto x_coord = domain.geom_traits().compute_x_3_object(); + auto y_coord = domain.geom_traits().compute_y_3_object(); + auto z_coord = domain.geom_traits().compute_z_3_object(); + auto point = domain.geom_traits().construct_point_3_object(); + + if((val_0 <= isovalue) == (val_1 <= isovalue)) + return false; + + const FT den = val_0 - val_1; + const FT u = is_zero(den) ? 0.5 : (val_0 - isovalue) / den; + p = point((FT(1) - u) * x_coord(p_0) + u * x_coord(p_1), + (FT(1) - u) * y_coord(p_0) + u * y_coord(p_1), + (FT(1) - u) * z_coord(p_0) + u * z_coord(p_1)); + + return true; + } +}; + +/* + * \ingroup IS_Domain_helpers_grp + * + * \cgalModels{EdgeIntersectionOracle_3} + * + * \brief The class `Ray_marching_edge_intersection` uses ray marching to find the intersection point + * between an edge and the isosurface. + * + * This class is suitable when the values stem from a signed distance function. + */ +struct Ray_marching_edge_intersection +{ + template + bool operator()(const typename Domain::Edge_descriptor& e, + const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + typename Domain::Geom_traits::Point_3& p) const + { + // @todo this is for the case where we know domain.value is an SDF + // then we can do better than a dichotomy + // Take code from the AW3 sharp branch + CGAL_assertion(false); + return false; + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_EDGE_INTERSECTION_ORACLES_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Cell_type.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Cell_type.h index ffd37f03bce..06ed30b7411 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Cell_type.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Cell_type.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -9,8 +9,8 @@ // // Author(s) : Julian Stahl -#ifndef CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE -#define CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE +#ifndef CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H +#define CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H #include @@ -19,16 +19,16 @@ namespace CGAL { namespace Isosurfacing { -// Was supposed to check if an algorithm can handle a specific domain. Not used right now. @fixme remove +// Was supposed to check if an algorithm can handle a specific domain. Not used right now. using Cell_type = std::size_t; -static constexpr Cell_type ANY_CELL = (std::numeric_limits::max)(); +static constexpr Cell_type ANY_CELL = (std::numeric_limits::max)(); -static constexpr Cell_type POLYHEDRAL_CELL = (Cell_type(1) << 0); -static constexpr Cell_type TETRAHEDRAL_CELL = (Cell_type(1) << 1); -static constexpr Cell_type CUBICAL_CELL = (Cell_type(1) << 2); +static constexpr Cell_type POLYHEDRAL_CELL = (std::size_t(1) << 0); +static constexpr Cell_type TETRAHEDRAL_CELL = (std::size_t(1) << 1); +static constexpr Cell_type CUBICAL_CELL = (std::size_t(1) << 2); } // namespace Isosurfacing } // namespace CGAL -#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE +#endif // CGAL_ISOSURFACING_3_INTERNAL_DOMAIN_CELL_TYPE_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h deleted file mode 100644 index 50ae0fc6ba9..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_function.h +++ /dev/null @@ -1,52 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H -#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H - -#include - -#include - -namespace CGAL { -namespace Isosurfacing { -namespace internal { - -template -class Explicit_Cartesian_grid_function -{ -public: - using Grid = Grid_; - using Geom_traits = typename Grid::Geom_traits; - using FT = typename Geom_traits::FT; - - using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor; - -public: - Explicit_Cartesian_grid_function(const Grid& grid) - : m_grid{grid} - { } - - // gets the value at vertex `v` - FT operator()(const Vertex_descriptor& v) const - { - return m_grid.value(v[0], v[1], v[2]); - } - -private: - const Grid& m_grid; -}; - -} // namespace internal -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_FUNCTION_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h deleted file mode 100644 index bbc090a2946..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Explicit_Cartesian_grid_geometry_3.h +++ /dev/null @@ -1,47 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H -#define CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H - -#include - -#include - -namespace CGAL { -namespace Isosurfacing { -namespace internal { - -template -class Explicit_Cartesian_grid_geometry_3 -{ - using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor; - -public: - Explicit_Cartesian_grid_geometry_3(const Grid& grid) - : m_grid{grid} - { } - - // gets the position of vertex `v` - decltype(auto) /*Point_3*/ operator()(const Vertex_descriptor& v) const - { - return m_grid.point(v[0], v[1], v[2]); - } - -private: - const Grid& m_grid; -}; - -} // namespace internal -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_INTERNAL_EXPLICIT_CARTESIAN_GRID_GEOMETRY_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h deleted file mode 100644 index bc8a61d8ca2..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_Cartesian_grid_geometry_3.h +++ /dev/null @@ -1,67 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H -#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H - -#include - -#include - -namespace CGAL { -namespace Isosurfacing { -namespace internal { - -// Describes the geometry of a regular Cartesian grid. -// Positions are not stored but calculated on-the-fly from an offset and grid spacing. -template -class Implicit_Cartesian_grid_geometry_3 -{ -public: - using Geom_traits = GeomTraits; - using Point_3 = typename Geom_traits::Point_3; - using Vector_3 = typename Geom_traits::Vector_3; - - using Vertex_descriptor = typename Grid_topology_3::Vertex_descriptor; - -private: - const Vector_3 m_offset; - const Vector_3 m_spacing; - Geom_traits m_gt; - -public: - Implicit_Cartesian_grid_geometry_3(const Vector_3& offset, - const Vector_3& spacing, - const Geom_traits& gt = Geom_traits()) - : m_offset{offset}, - m_spacing{spacing}, - m_gt{gt} - { } - - // gets the position of vertex v - Point_3 operator()(const Vertex_descriptor& v) const - { - typename Geom_traits::Compute_x_3 x_coord = m_gt.compute_x_3_object(); - typename Geom_traits::Compute_y_3 y_coord = m_gt.compute_y_3_object(); - typename Geom_traits::Compute_z_3 z_coord = m_gt.compute_z_3_object(); - typename Geom_traits::Construct_point_3 point = m_gt.construct_point_3_object(); - - return point(v[0] * x_coord(m_spacing) + x_coord(m_offset), - v[1] * y_coord(m_spacing) + y_coord(m_offset), - v[2] * z_coord(m_spacing) + z_coord(m_offset)); - } -}; - -} // namespace internal -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_CARTESIAN_GRID_GEOMETRY_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h deleted file mode 100644 index 7493e168a69..00000000000 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h +++ /dev/null @@ -1,53 +0,0 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// Author(s) : Julian Stahl - -#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H -#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H - -#include - -namespace CGAL { -namespace Isosurfacing { -namespace internal { - -// Wrapper for an implicit function that can only be evaluated at a position and not at a vertex. -// Evaluates the geometry to get the vertex position and passes that to the function. -template -class Implicit_function_with_geometry -{ - using Point_function = PointFunction; - -public: - // creates a function that uses the geometry to evaluate the function at vertex positions. - Implicit_function_with_geometry(const Geometry& geom, - const Point_function& func) - : m_geom(geom), - m_func(func) - { } - - // gets the value of the function at vertex `v` - template - decltype(auto) /*FT*/ operator()(const VertexDescriptor& v) const - { - return m_func.operator()(m_geom.operator()(v)); - } - -private: - const Geometry m_geom; - const Point_function m_func; -}; - -} // namespace internal -} // namespace Isosurfacing -} // namespace CGAL - -#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_FUNCTION_WITH_GEOMETRY_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h index 5b86dfc9b35..9c91eb2ab47 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Isosurfacing_domain_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,140 +8,153 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Julian Stahl +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H #define CGAL_ISOSURFACING_3_INTERNAL_ISOSURFACING_DOMAIN_3_H #include -#include +#include +#include + +#include namespace CGAL { namespace Isosurfacing { namespace internal { -// A wrapper class to puzzle a domain together from different combinations of topology, -// geometry, function, and gradient. -template +// This class is pretty much just the concatenation of the following classes: +// - Partition: Space partitioning data structure, e.g. Cartesian grid, octree, ... +// - Values: values over the 3D space +// - Gradients: gradients over the 3D space +// - Oracle: edge-isosurface intersection computation +template class Isosurfacing_domain_3 { public: - using Geom_traits = GeomTraits; + using Edge_intersection_oracle = IntersectionOracle; + using Geom_traits = typename Partition::Geom_traits; using FT = typename Geom_traits::FT; using Point_3 = typename Geom_traits::Point_3; using Vector_3 = typename Geom_traits::Vector_3; - using Topology = Topology_; - using Vertex_descriptor = typename Topology_::Vertex_descriptor; - using Edge_descriptor = typename Topology_::Edge_descriptor; - using Cell_descriptor = typename Topology_::Cell_descriptor; + using PT = CGAL::Isosurfacing::partition_traits; - using Vertices_incident_to_edge = typename Topology_::Vertices_incident_to_edge; - using Cells_incident_to_edge = typename Topology_::Cells_incident_to_edge; - using Cell_vertices = typename Topology_::Cell_vertices; - using Cell_edges = typename Topology_::Cell_edges; + using Vertex_descriptor = typename PT::Vertex_descriptor; + using Edge_descriptor = typename PT::Edge_descriptor; + using Cell_descriptor = typename PT::Cell_descriptor; - using Geometry = Geometry_; - using Function = Function_; - using Gradient = Gradient_; + using Vertices_incident_to_edge = typename PT::Vertices_incident_to_edge; + using Cells_incident_to_edge = typename PT::Cells_incident_to_edge; + using Cell_vertices = typename PT::Cell_vertices; + using Cell_edges = typename PT::Cell_edges; - static constexpr Cell_type CELL_TYPE = Topology_::CELL_TYPE; - static constexpr std::size_t VERTICES_PER_CELL = Topology_::VERTICES_PER_CELL; - static constexpr std::size_t EDGES_PER_CELL = Topology_::EDGES_PER_CELL; + static constexpr Cell_type CELL_TYPE = PT::CELL_TYPE; + static constexpr std::size_t VERTICES_PER_CELL = PT::VERTICES_PER_CELL; + static constexpr std::size_t EDGES_PER_CELL = PT::EDGES_PER_CELL; private: - const Topology m_topo; - const Geometry m_geom; - const Function m_func; - const Gradient m_grad; - const Geom_traits m_gt; + const Partition& m_partition; + const ValueField& m_values; + const GradientField& m_gradients; + const IntersectionOracle m_intersection_oracle; public: - // creates a base domain from topology, geometry, implicit function, gradient, and geometric traits - Isosurfacing_domain_3(const Topology& topo, - const Geometry& geom, - const Function& func, - const Gradient& grad, - const Geom_traits& gt) - : m_topo{topo}, - m_geom{geom}, - m_func{func}, - m_grad{grad}, - m_gt(gt) + Isosurfacing_domain_3(const Partition& partition, + const ValueField& values, + const GradientField& gradients, + const IntersectionOracle& intersection_oracle = IntersectionOracle()) + : m_partition{partition}, + m_values{values}, + m_gradients{gradients}, + m_intersection_oracle{intersection_oracle} { } const Geom_traits& geom_traits() const { - return m_gt; + return m_partition.geom_traits(); + } + + const Edge_intersection_oracle& intersection_oracle() const + { + return m_intersection_oracle; } public: + // The following functions are dispatching to the partition_traits' static functions. + // gets the position of vertex `v` decltype(auto) /*Point_3*/ point(const Vertex_descriptor& v) const { - return m_geom.operator()(v); - } - - // gets the value of the function at vertex `v` - decltype(auto) /*FT*/ value(const Vertex_descriptor& v) const - { - return m_func.operator()(v); - } - - // gets the gradient at vertex `v` - decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const - { - return m_grad.operator()(p); + return PT::point(v, m_partition); } // gets a container with the two vertices incident to the edge `e` decltype(auto) /*Vertices_incident_to_edge*/ incident_vertices(const Edge_descriptor& e) const { - return m_topo.incident_vertices(e); + return PT::incident_vertices(e, m_partition); } // gets a container with all cells incident to the edge `e` decltype(auto) /*Cells_incident_to_edge*/ incident_cells(const Edge_descriptor& e) const { - return m_topo.incident_cells(e); + return PT::incident_cells(e, m_partition); } // gets a container with all vertices of the cell `c` decltype(auto) /*Cell_vertices*/ cell_vertices(const Cell_descriptor& c) const { - return m_topo.cell_vertices(c); + return PT::cell_vertices(c, m_partition); } // gets a container with all edges of the cell `c` decltype(auto) /*Cell_edges*/ cell_edges(const Cell_descriptor& c) const { - return m_topo.cell_edges(c); + return PT::cell_edges(c, m_partition); } // iterates over all vertices `v`, calling `f(v)` on each of them template void for_each_vertex(Functor& f) const { - m_topo.for_each_vertex(f, ConcurrencyTag{}); + PT::for_each_vertex(f, m_partition, ConcurrencyTag{}); } // iterates over all edges `e`, calling `f(e)` on each of them template void for_each_edge(Functor& f) const { - m_topo.for_each_edge(f, ConcurrencyTag{}); + PT::for_each_edge(f, m_partition, ConcurrencyTag{}); } // iterates over all cells `c`, calling `f(c)` on each of them template void for_each_cell(Functor& f) const { - m_topo.for_each_cell(f, ConcurrencyTag{}); + PT::for_each_cell(f, m_partition, ConcurrencyTag{}); + } + + // gets the value of the function at vertex `v` + decltype(auto) /*FT*/ value(const Vertex_descriptor& v) const + { + return m_values(v); + } + + // gets the value of the function at point `p` + decltype(auto) /*FT*/ value(const Point_3& p) const + { + return m_values(p); + } + + // gets the gradient at point `p` + decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const + { + return m_gradients(p); } }; diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h index 51bb87365f4..67813322e41 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h index 90e82c0a4d9..7128a9b269b 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_wrapper.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_wrapper.h index 6b81e494f37..373804672e3 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_wrapper.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_wrapper.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -30,6 +31,7 @@ namespace internal { template class Octree_wrapper + : public Octree_topology > > { /* * Naming convention from "A parallel dual marching cubes approach to quad diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h index 9fdae159c15..e9c93f3953a 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/dual_contouring_functors.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -9,13 +9,15 @@ // // Author(s) : Daniel Zint // Julian Stahl +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_FUNCTORS_H #define CGAL_ISOSURFACING_3_INTERNAL_DUAL_CONTOURING_FUNCTORS_H #include -#include +#include +#include #include #ifdef CGAL_EIGEN3_ENABLED @@ -25,6 +27,14 @@ #include #endif +#ifdef CGAL_LINKED_WITH_TBB +#if TBB_INTERFACE_VERSION < 12010 && !defined(TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS) +#define CGAL_HAS_DEFINED_TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS +#define TBB_PREVIEW_CONCURRENT_ORDERED_CONTAINERS 1 +#endif +#include +#endif +#include #include #include #include @@ -32,406 +42,616 @@ namespace CGAL { 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. - * - * \tparam constrain_to_cell clamp vertex position to the bounding box of the cell - */ -template -class QEM_SVD -{ -public: - /* - * \brief computes the vertex position for a point in Dual Contouring. - * - * \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`. - * - * \param domain the domain providing input data and its topology - * \param isovalue value of the isosurface - * \param cell the cell within the domain for which the vertex position is computed - * \param point the point position of the vertex that belongs to that cell - * - * \return `true` if the voxel intersects the isosurface, `false` otherwise - */ - template - bool position(const Domain& domain, - const typename Domain::Geom_traits::FT isovalue, - const typename Domain::Cell_descriptor& cell, - typename Domain::Geom_traits::Point_3& p) const - { - using Geom_traits = typename Domain::Geom_traits; - using FT = typename Geom_traits::FT; - using Point_3 = typename Geom_traits::Point_3; - using Vector_3 = typename Geom_traits::Vector_3; - - using Eigen_vector_3 = Eigen_vector; - using Eigen_vector_x = Eigen_vector; - using Eigen_matrix_3 = Eigen_matrix; - using Eigen_matrix_x = Eigen_matrix; - - using Vertex_descriptor = typename Domain::Vertex_descriptor; - - auto x_coord = domain.geom_traits().compute_x_3_object(); - auto y_coord = domain.geom_traits().compute_y_3_object(); - auto z_coord = domain.geom_traits().compute_z_3_object(); - auto point = domain.geom_traits().construct_point_3_object(); - - typename Domain::Cell_vertices vertices = domain.cell_vertices(cell); - const std::size_t cn = vertices.size(); - - // compute edge intersections - std::vector edge_intersections; - std::vector edge_intersection_normals; - - for(const auto& edge : domain.cell_edges(cell)) - { - const auto& edge_vertices = domain.incident_vertices(edge); - const Vertex_descriptor& v0 = edge_vertices[0]; - const Vertex_descriptor& v1 = edge_vertices[1]; - - const FT val0 = domain.value(v0); - const FT val1 = domain.value(v1); - - const Point_3& p0 = domain.point(v0); - const Point_3& p1 = domain.point(v1); - - if((val0 <= isovalue) != (val1 <= isovalue)) - { - // current edge is intersected by the isosurface - const FT u = (val0 - isovalue) / (val0 - val1); - const Point_3 p_lerp = point((FT(1) - u) * x_coord(p0) + u * x_coord(p1), - (FT(1) - u) * y_coord(p0) + u * y_coord(p1), - (FT(1) - u) * z_coord(p0) + u * z_coord(p1)); - edge_intersections.push_back(p_lerp); - edge_intersection_normals.push_back(domain.gradient(p_lerp)); - } - } - - if(edge_intersections.empty()) - return false; - - FT x_min, y_min, z_min, x_max, y_max, z_max; - FT x(0), y(0), z(0); - - for(const auto& v : vertices) - { - const Point_3 cp = domain.point(v); - x += x_coord(cp); - y += y_coord(cp); - z += z_coord(cp); - - if constexpr(constrain_to_cell) - { - x_min = (std::min)(x_min, x_coord(cp)); - y_min = (std::min)(y_min, y_coord(cp)); - z_min = (std::min)(z_min, z_coord(cp)); - - x_max = (std::max)(x_max, x_coord(cp)); - y_max = (std::max)(y_max, y_coord(cp)); - z_max = (std::max)(z_max, z_coord(cp)); - } - } - - p = point(x / cn, y / cn, z / cn); - - // SVD QEM - Eigen_matrix_3 A; - A.setZero(); - Eigen_vector_3 rhs; - rhs.setZero(); - for(std::size_t i=0; i etc., - // so the double conversion is not implicit - Eigen_matrix_3 A_k = typename Eigen_matrix_3::EigenType(n_k * n_k.transpose()); - Eigen_vector_3 b_k; - b_k = d_k * n_k; - A += A_k; - rhs += b_k; - } - - Eigen::JacobiSVD 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_vector_3 x_hat; - x_hat << x_coord(p), y_coord(p), z_coord(p); - - // Lindstrom formula for QEM new position for singular matrices - Eigen_vector_x v_svd; - v_svd = x_hat + svd.solve(rhs - A * x_hat); - - // bounding box - if constexpr(constrain_to_cell) - { - v_svd[0] = std::clamp(v_svd[0], x_min, x_max); - v_svd[1] = std::clamp(v_svd[1], y_min, y_max); - v_svd[2] = std::clamp(v_svd[2], z_min, z_max); - } - - p = point(v_svd[0], v_svd[1], v_svd[2]); - - return true; - } -}; - -/** - * \brief returns the cell's center. - */ -class Cell_center -{ -public: - /* - * \brief computes the vertex position for a point in Dual Contouring. - * - * \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`. - * - * \param domain the domain providing input data and its topology - * \param isovalue value of the isosurface - * \param cell the cell within the domain for which the vertex position is computed - * \param point the point position of the vertex that belongs to that cell - * - * \return `true` if the voxel intersects the isosurface, `false` otherwise - */ - template - bool position(const Domain& domain, - const typename Domain::Geom_traits::FT isovalue, - const typename Domain::Cell_descriptor& vh, - typename Domain::Geom_traits::Point_3& p) const - { - using FT = typename Domain::Geom_traits::FT; - using Point_3 = typename Domain::Geom_traits::Vector_3; - - using Vertex_descriptor = typename Domain::Vertex_descriptor; - - auto x_coord = domain.geom_traits().compute_x_3_object(); - auto y_coord = domain.geom_traits().compute_y_3_object(); - auto z_coord = domain.geom_traits().compute_z_3_object(); - auto point = domain.geom_traits().construct_point_3_object(); - - typename Domain::Cell_vertices vertices = domain.cell_vertices(vh); - const std::size_t cn = vertices.size(); - - bool all_smaller = true; - bool all_greater = true; - for(const auto& v : vertices) - { - const bool b = (domain.value(v) <= isovalue); - all_smaller = all_smaller && b; - all_greater = all_greater && !b; - } - - if(all_smaller || all_greater) - return false; - - FT x(0), y(0), z(0); - for(const auto& v : vertices) - { - const Point_3 cp = domain.point(v); - x += x_coord(cp); - y += y_coord(cp); - z += z_coord(cp); - } - - // set point to cell center - p = point(x / cn, y / cn, z / cn); - - return true; - } -}; - -/* - * \brief computes the centroid of all cell edge intersections with the isosurface. - */ -class Centroid_of_edge_intersections -{ -public: - /* - * \brief computes the vertex position for a point in Dual Contouring. - * - * \tparam Domain must be a model of `IsosurfacingDomainWithGradient_3`. - * - * \param domain the domain providing input data and its topology - * \param isovalue value of the isosurface - * \param cell the cell within the domain for which the vertex position is computed - * \param point the point position of the vertex that belongs to that cell - * - * \return `true` if the voxel intersects the isosurface, `false` otherwise - */ - template - bool position(const Domain& domain, - const typename Domain::Geom_traits::FT isovalue, - const typename Domain::Cell_descriptor& cell, - typename Domain::Geom_traits::Point_3& p) const - { - using FT = typename Domain::Geom_traits::FT; - using Point_3 = typename Domain::Geom_traits::Point_3; - using Vector_3 = typename Domain::Geom_traits::Vector_3; - - using Vertex_descriptor = typename Domain::Vertex_descriptor; - using Edge_descriptor = typename Domain::Edge_descriptor; - - auto x_coord = domain.geom_traits().compute_x_3_object(); - auto y_coord = domain.geom_traits().compute_y_3_object(); - auto z_coord = domain.geom_traits().compute_z_3_object(); - auto point = domain.geom_traits().construct_point_3_object(); - - // compute edge intersections - std::vector edge_intersections; - - for(const Edge_descriptor& edge : domain.cell_edges(cell)) - { - const auto& edge_vertices = domain.incident_vertices(edge); - const Vertex_descriptor& v0 = edge_vertices[0]; - const Vertex_descriptor& v1 = edge_vertices[1]; - - const FT val_0 = domain.value(v0); - const FT val_1 = domain.value(v1); - - const Point_3& p0 = domain.point(v0); - const Point_3& p1 = domain.point(v1); - - if((val_0 <= isovalue) != (val_1 <= isovalue)) - { - // current edge is intersected by the isosurface - const FT u = (val_0 - isovalue) / (val_0 - val_1); - const Point_3 p_lerp = point((FT(1) - u) * x_coord(p0) + u * x_coord(p1), - (FT(1) - u) * y_coord(p0) + u * y_coord(p1), - (FT(1) - u) * z_coord(p0) + u * z_coord(p1)); - edge_intersections.push_back(p_lerp); - } - } - - const std::size_t en = edge_intersections.size(); - if(en == 0) - return false; - - FT x(0), y(0), z(0); - for(const Point_3& p : edge_intersections) - { - x += x_coord(p); - y += y_coord(p); - z += z_coord(p); - } - - p = point(x / en, y / en, z / en); - - return true; - } -}; - -} // namespace Positioning template -class Dual_contouring_vertex_positioning + typename EdgeToPointIDMap, + typename PointRange, + typename GradientRange> +bool cell_position_QEM(const typename Domain::Cell_descriptor& c, + const Domain& domain, + const bool constrain_to_cell, + const EdgeToPointIDMap& edge_to_point_id, + const PointRange& edge_points, + const GradientRange& edge_gradients, + typename Domain::Geom_traits::Point_3& p) +{ + using Geom_traits = typename Domain::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; + + using Eigen_vector_3 = Eigen_vector; + using Eigen_matrix_3 = Eigen_matrix; + using Eigen_vector_x = Eigen_vector; + using Eigen_matrix_x = Eigen_matrix; + + typename Geom_traits::Compute_x_3 x_coord = domain.geom_traits().compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = domain.geom_traits().compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = domain.geom_traits().compute_z_3_object(); + typename Geom_traits::Construct_point_3 point = domain.geom_traits().construct_point_3_object(); + + // compute edge intersections + std::vector cell_edge_intersections; + std::vector cell_edge_intersection_normals; + + for(const auto& edge : domain.cell_edges(c)) + { + const auto it = edge_to_point_id.find(edge); + if(it == edge_to_point_id.end()) + continue; + + cell_edge_intersections.push_back(edge_points[it->second]); + cell_edge_intersection_normals.push_back(edge_gradients[it->second]); + } + + const std::size_t en = cell_edge_intersections.size(); + if(en == 0) + return false; + +#ifdef CGAL_ISOSURFACING_3_DC_FUNCTORS_DEBUG + std::cout << "Points and normals: " << std::endl; + for(std::size_t i=0; i::max(); // @todo domain.bbox() + x_max = y_max = z_max = std::numeric_limits::lowest(); + FT x(0), y(0), z(0); + + typename Domain::Cell_vertices vertices = domain.cell_vertices(c); + for(const auto& v : vertices) + { + const Point_3& cp = domain.point(v); + if(constrain_to_cell) + { + x_min = (std::min)(x_min, x_coord(cp)); + y_min = (std::min)(y_min, y_coord(cp)); + z_min = (std::min)(z_min, z_coord(cp)); + + x_max = (std::max)(x_max, x_coord(cp)); + y_max = (std::max)(y_max, y_coord(cp)); + z_max = (std::max)(z_max, z_coord(cp)); + } + } + + for(const auto& ep : cell_edge_intersections) + { + x += x_coord(ep); + y += y_coord(ep); + z += z_coord(ep); + } + + Point_3 com = point(x / en, y / en, z / en); + +#ifdef CGAL_ISOSURFACING_3_DC_FUNCTORS_DEBUG + std::cout << "cell: " << x_min << " " << y_min << " " << z_min << " " << x_max << " " << y_max << " " << z_max << std::endl; + std::cout << "COM " << com << std::endl; +#endif + + // SVD QEM + Eigen_matrix_3 A; + A.setZero(); + Eigen_vector_3 rhs; + rhs.setZero(); + for(std::size_t i=0; i etc., + // so the double conversion is not implicit + Eigen_matrix_3 A_k = typename Eigen_matrix_3::EigenType(n_k * n_k.transpose()); + Eigen_vector_3 b_k; + b_k = d_k * n_k; + A += A_k; + rhs += b_k; + } + + Eigen::JacobiSVD 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); + + Eigen_vector_3 x_hat; + x_hat << x_coord(com), y_coord(com), z_coord(com); + + // Lindstrom formula for QEM new position for singular matrices + Eigen_vector_x v_svd; + v_svd = x_hat + svd.solve(rhs - A * x_hat); + + + // bounding box + if(constrain_to_cell) + { + v_svd[0] = std::clamp(v_svd[0], x_min, x_max); + v_svd[1] = std::clamp(v_svd[1], y_min, y_max); + v_svd[2] = std::clamp(v_svd[2], z_min, z_max); + } + + p = point(v_svd[0], v_svd[1], v_svd[2]); + +#ifdef CGAL_ISOSURFACING_3_DC_FUNCTORS_DEBUG + std::cout << "CGAL QEM POINT: " << v_svd[0] << " " << v_svd[1] << " " << v_svd[2] << std::endl; + std::cout << "CGAL clamped QEM POINT: " << p[0] << " " << p[1] << " " << p[2] << std::endl; + std::cout << "--" << std::endl; +#endif + + return true; +} + +template +void generate_face(const typename Domain::Edge_descriptor& e, + const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + const bool do_not_triangulate_faces, + const EdgeToPointIDMap& edge_to_point_id, + const CellToPointIDMap& cell_to_point_id, + std::mutex& mutex, + PolygonRange& polygons) { -private: using FT = typename Domain::Geom_traits::FT; - using Point_3 = typename Domain::Geom_traits::Point_3; using Cell_descriptor = typename Domain::Cell_descriptor; -public: - Dual_contouring_vertex_positioning(const Domain& domain, - const FT isovalue, - const Positioning& positioning) - : domain(domain), - isovalue(isovalue), - positioning(positioning), - points_counter(0) - { } + const auto& vertices = domain.incident_vertices(e); - void operator()(const Cell_descriptor& v) + // @todo this check could be avoided for QEM: active edges are in `edge_to_point_id` + const FT val_0 = domain.value(vertices[0]); + const FT val_1 = domain.value(vertices[1]); + if((val_0 <= isovalue) == (val_1 <= isovalue)) + return; + + std::vector vertex_ids; + + const auto& icells = domain.incident_cells(e); + for(const Cell_descriptor& c : icells) { - // compute dc-vertices - Point_3 p; // @fixme initialize? - if(positioning.position(domain, isovalue, v, p)) - { - std::lock_guard lock(mutex); - map_voxel_to_point[v] = p; - map_voxel_to_point_id[v] = points_counter++; - } + auto it = cell_to_point_id.find(c); + if(it == cell_to_point_id.end()) + continue; + + vertex_ids.push_back(it->second); } - // private: // @todo - const Domain& domain; - const FT isovalue; - const Positioning& positioning; + if(vertex_ids.size() < 3) + return; - std::map map_voxel_to_point_id; // @todo unordered? - std::map map_voxel_to_point; - std::size_t points_counter; + if(val_0 > val_1) + std::reverse(vertex_ids.begin(), vertex_ids.end()); - std::mutex mutex; -}; + // @todo? filter degenerate faces? -template -class Dual_contouring_face_generation + if(do_not_triangulate_faces) + { + std::lock_guard lock(mutex); + polygons.push_back(vertex_ids); + } + else + { + auto it = edge_to_point_id.find(e); + if(it == edge_to_point_id.end()) + { + CGAL_assertion(false); + return; + } + const std::size_t ei = it->second; + + std::lock_guard lock(mutex); + for(std::size_t i=0; i +class Dual_contourer; + +template +class Dual_contourer { -private: using FT = typename Domain::Geom_traits::FT; + using Point_3 = typename Domain::Geom_traits::Point_3; + using Vector_3 = typename Domain::Geom_traits::Vector_3; + using Vertex_descriptor = typename Domain::Vertex_descriptor; using Edge_descriptor = typename Domain::Edge_descriptor; using Cell_descriptor = typename Domain::Cell_descriptor; + std::mutex m_mutex; + public: - Dual_contouring_face_generation(const Domain& domain, - const FT isovalue) - : domain(domain), - isovalue(isovalue) - { } - - void operator()(const Edge_descriptor& e) + template + void operator()(const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + PointRange& points, + PolygonRange& polygons, + const NamedParameters& np = parameters::default_values()) { - // save all faces - const auto& vertices = domain.incident_vertices(e); - const FT val_0 = domain.value(vertices[0]); - const FT val_1 = domain.value(vertices[1]); + using parameters::choose_parameter; + using parameters::get_parameter; - if(val_0 <= isovalue && val_1 > isovalue) + // Otherwise the `edge_to_point_id` map might be messed up + CGAL_precondition(points.empty()); + CGAL_precondition(polygons.empty()); + + const bool constrain_to_cell = choose_parameter(get_parameter(np, internal_np::constrain_to_cell), false); + + const bool do_not_triangulate_faces = + choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false); + + using Edge_to_point_ID_map = std::map; + using Cell_to_point_ID_map = std::map; + + Edge_to_point_ID_map edge_to_point_id; + Cell_to_point_ID_map cell_to_point_id; + + std::vector edge_points; + std::vector edge_gradients; + + // --------------------------------------------------------------------------------------------- + // construct the intersection of the surface at active edges + auto edge_positioner = [&](const Edge_descriptor& e) { - const auto& voxels = domain.incident_cells(e); + const auto& evs = domain.incident_vertices(e); + const Vertex_descriptor& v0 = evs[0]; + const Vertex_descriptor& v1 = evs[1]; + const Point_3& p0 = domain.point(v0); + const Point_3& p1 = domain.point(v1); + const FT val0 = domain.value(v0); + const FT val1 = domain.value(v1); - std::lock_guard lock(mutex); - faces[e].insert(faces[e].begin(), voxels.begin(), voxels.end()); - } - else if(val_1 <= isovalue && val_0 > isovalue) + Point_3 p; + bool res = domain.intersection_oracle()(p0, p1, val0, val1, domain, isovalue, p); + if(!res) + return; + + const Vector_3 g = domain.gradient(p); + + std::lock_guard lock(m_mutex); + edge_to_point_id[e] = edge_points.size(); + edge_points.push_back(p); + edge_gradients.push_back(g); + }; + domain.template for_each_edge(edge_positioner); + +#ifdef CGAL_ISOSURFACING_3_DC_FUNCTORS_DEBUG + auto x_coord = domain.geom_traits().compute_x_3_object(); + auto y_coord = domain.geom_traits().compute_y_3_object(); + auto z_coord = domain.geom_traits().compute_z_3_object(); + + std::ofstream out_active_edges("active_edges.polylines"); + for(const auto& ei : edge_to_point_id) { - const auto& voxels = domain.incident_cells(e); + const Edge_descriptor& e = ei.first; + const auto& evs = domain.incident_vertices(e); + const Vertex_descriptor& v0 = evs[0]; + const Vertex_descriptor& v1 = evs[1]; + const Point_3& p0 = domain.point(v0); + const Point_3& p1 = domain.point(v1); - std::lock_guard lock(mutex); - faces[e].insert(faces[e].begin(), voxels.rbegin(), voxels.rend()); + out_active_edges << "2 " << x_coord(p0) << " " << y_coord(p0) << " " << z_coord(p0) << " " << x_coord(p1) << " " << y_coord(p1) << " " << z_coord(p1) << std::endl; } + + std::ofstream out_edge_intersections("edge_intersections.polylines"); + for(const auto& ei : edge_to_point_id) + { + const Point_3& p = edge_points.at(ei.second); + const Vector_3& g = edge_gradients.at(ei.second); + + out_edge_intersections << "2 " << x_coord(p) << " " << y_coord(p) << " " << z_coord(p) << " " << x_coord(p) + x_coord(g) << " " << y_coord(p) + y_coord(g) << " " << z_coord(p) + z_coord(g) << std::endl; + } + +#endif + + if(!do_not_triangulate_faces) + points.insert(points.end(), edge_points.begin(), edge_points.end()); + + // --------------------------------------------------------------------------------------------- + // create a vertex for each cell that has at least one active edge + auto cell_positioner = [&](const Cell_descriptor& c) + { + Point_3 p; + if(cell_position_QEM(c, domain, constrain_to_cell, edge_to_point_id, + edge_points, edge_gradients, p)) + { + std::lock_guard lock(m_mutex); // @todo useless if sequential + cell_to_point_id[c] = points.size(); + points.push_back(p); + } + }; + domain.template for_each_cell(cell_positioner); + + // --------------------------------------------------------------------------------------------- + // connect vertices around edges to form faces + auto face_generator = [&](const Edge_descriptor& e) + { + generate_face(e, domain, isovalue, do_not_triangulate_faces, + edge_to_point_id, cell_to_point_id, m_mutex, polygons); + }; + domain.template for_each_edge(face_generator); + +#ifdef CGAL_ISOSURFACING_3_DC_FUNCTORS_DEBUG + std::cout << points.size() << " points" << std::endl; + std::cout << polygons.size() << " polygons" << std::endl; +#endif } +}; - // private: // @todo - std::map > faces; +template +class Dual_contourer +{ + using Geom_traits = typename Domain::Geom_traits; - const Domain& domain; - const FT isovalue; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; - std::mutex mutex; + using Vertex_descriptor = typename Domain::Vertex_descriptor; + using Edge_descriptor = typename Domain::Edge_descriptor; + using Cell_descriptor = typename Domain::Cell_descriptor; + + std::mutex m_mutex; + +public: + template + void operator()(const Domain& domain, + const typename Geom_traits::FT isovalue, + PointRange& points, + PolygonRange& polygons, + const NamedParameters& np = parameters::default_values()) + { + using parameters::choose_parameter; + using parameters::get_parameter; + + // Otherwise the `edge_to_point_id` map might be messed up + CGAL_precondition(points.empty()); + CGAL_precondition(polygons.empty()); + + bool do_not_triangulate_faces = + choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false); + + using Edge_to_point_ID_map = std::map; + using Cell_to_point_ID_map = std::map; + + Edge_to_point_ID_map edge_to_point_id; + Cell_to_point_ID_map cell_to_point_id; + + std::vector edge_points; + + // --------------------------------------------------------------------------------------------- + auto edge_positioner = [&](const Edge_descriptor& e) + { + const auto& evs = domain.incident_vertices(e); + const Vertex_descriptor& v0 = evs[0]; + const Vertex_descriptor& v1 = evs[1]; + const Point_3& p0 = domain.point(v0); + const Point_3& p1 = domain.point(v1); + const FT val0 = domain.value(v0); + const FT val1 = domain.value(v1); + + Point_3 p; + bool res = domain.intersection_oracle()(p0, p1, val0, val1, domain, isovalue, p); + if(!res) + return; + + std::lock_guard lock(m_mutex); + edge_to_point_id[e] = edge_points.size(); + edge_points.push_back(p); + }; + domain.template for_each_edge(edge_positioner); + + if(!do_not_triangulate_faces) + points.insert(points.end(), edge_points.begin(), edge_points.end()); + + // --------------------------------------------------------------------------------------------- + auto cell_positioner = [&](const Cell_descriptor& c) + { + typename Geom_traits::Compute_x_3 x_coord = domain.geom_traits().compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = domain.geom_traits().compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = domain.geom_traits().compute_z_3_object(); + typename Geom_traits::Construct_point_3 point = domain.geom_traits().construct_point_3_object(); + + // compute edge intersections + std::vector edge_intersections; + for(const Edge_descriptor& e : domain.cell_edges(c)) + { + const auto it = edge_to_point_id.find(e); + if(it == edge_to_point_id.end()) + continue; + + edge_intersections.push_back(edge_points[it->second]); // @todo could avoid copying + } + + const std::size_t en = edge_intersections.size(); + if(en == 0) + return; + + FT x = 0, y = 0, z = 0; + for(const Point_3& p : edge_intersections) + { + x += x_coord(p); + y += y_coord(p); + z += z_coord(p); + } + + const Point_3 p = point(x / en, y / en, z / en); + + std::lock_guard lock(m_mutex); + cell_to_point_id[c] = points.size(); + points.push_back(p); + }; + domain.template for_each_cell(cell_positioner); + + // --------------------------------------------------------------------------------------------- + auto face_generator = [&](const Edge_descriptor& e) + { + generate_face(e, domain, isovalue, do_not_triangulate_faces, + edge_to_point_id, cell_to_point_id, m_mutex, polygons); + }; + domain.template for_each_edge(face_generator); + } +}; + +template +class Dual_contourer +{ + using FT = typename Domain::Geom_traits::FT; + using Point_3 = typename Domain::Geom_traits::Point_3; + using Vector_3 = typename Domain::Geom_traits::Vector_3; + + using Vertex_descriptor = typename Domain::Vertex_descriptor; + using Edge_descriptor = typename Domain::Edge_descriptor; + using Cell_descriptor = typename Domain::Cell_descriptor; + + std::mutex m_mutex; + +public: + template + void operator()(const Domain& domain, + const typename Domain::Geom_traits::FT isovalue, + PointRange& points, + PolygonRange& polygons, + const NamedParameters& np = parameters::default_values()) + { + using parameters::choose_parameter; + using parameters::get_parameter; + + // Otherwise the `edge_to_point_id` map might be messed up + CGAL_precondition(points.empty()); + CGAL_precondition(polygons.empty()); + + bool do_not_triangulate_faces = + choose_parameter(get_parameter(np, internal_np::do_not_triangulate_faces), false); + + using Edge_to_point_ID_map = std::map; + using Cell_to_point_ID_map = std::map; + + Edge_to_point_ID_map edge_to_point_id; + Cell_to_point_ID_map cell_to_point_id; + + std::vector edge_points; + + // --------------------------------------------------------------------------------------------- + auto edge_positioner = [&](const Edge_descriptor& e) + { + using FT = typename Domain::Geom_traits::FT; + using Point_3 = typename Domain::Geom_traits::Point_3; + + using Vertex_descriptor = typename Domain::Vertex_descriptor; + + auto x_coord = domain.geom_traits().compute_x_3_object(); + auto y_coord = domain.geom_traits().compute_y_3_object(); + auto z_coord = domain.geom_traits().compute_z_3_object(); + auto point = domain.geom_traits().construct_point_3_object(); + + const auto& evs = domain.incident_vertices(e); + const Vertex_descriptor& v0 = evs[0]; + const Vertex_descriptor& v1 = evs[1]; + + const FT val_0 = domain.value(v0); + const FT val_1 = domain.value(v1); + if((val_0 <= isovalue) == (val_1 <= isovalue)) + return; + + Point_3 p = point((x_coord(domain.point(v0)) + x_coord(domain.point(v1))) / FT(2), + (y_coord(domain.point(v0)) + y_coord(domain.point(v1))) / FT(2), + (z_coord(domain.point(v0)) + z_coord(domain.point(v1))) / FT(2)); + + std::lock_guard lock(m_mutex); + edge_to_point_id[e] = edge_points.size(); + edge_points.push_back(p); + }; + domain.template for_each_edge(edge_positioner); + + if(!do_not_triangulate_faces) + points.insert(points.end(), edge_points.begin(), edge_points.end()); + + // --------------------------------------------------------------------------------------------- + auto cell_positioner = [&](const Cell_descriptor& c) + { + auto x_coord = domain.geom_traits().compute_x_3_object(); + auto y_coord = domain.geom_traits().compute_y_3_object(); + auto z_coord = domain.geom_traits().compute_z_3_object(); + auto point = domain.geom_traits().construct_point_3_object(); + + typename Domain::Cell_vertices vertices = domain.cell_vertices(c); + const std::size_t cn = vertices.size(); + + bool all_smaller = true; + bool all_greater = true; + for(const auto& v : vertices) + { + const bool b = (domain.value(v) <= isovalue); + all_smaller = all_smaller && b; + all_greater = all_greater && !b; + } + + if(all_smaller || all_greater) + return; + + FT x(0), y(0), z(0); + for(const auto& v : vertices) + { + const Point_3& cp = domain.point(v); + x += x_coord(cp); + y += y_coord(cp); + z += z_coord(cp); + } + + // set point to cell center + Point_3 p = point(x / cn, y / cn, z / cn); + + std::lock_guard lock(m_mutex); + cell_to_point_id[c] = points.size(); + points.push_back(p); + }; + domain.template for_each_cell(cell_positioner); + + // --------------------------------------------------------------------------------------------- + auto face_generator = [&](const Edge_descriptor& e) + { + generate_face(e, domain, isovalue, do_not_triangulate_faces, + edge_to_point_id, cell_to_point_id, m_mutex, polygons); + }; + domain.template for_each_edge(face_generator); + } }; } // namespace internal diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h new file mode 100644 index 00000000000..ea7b7446b3b --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h @@ -0,0 +1,182 @@ +// Copyright (c) 2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H +#define CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H + +#include +#include + +namespace CGAL { +namespace Isosurfacing { +namespace Shapes { + +// Shapes are defined at isovalue 0 + +// c is the center +// r the radius +template +typename K::FT +sphere(const typename K::Point_3& c, + const typename K::FT r, + const typename K::Point_3& q) +{ + return CGAL::approximate_sqrt(CGAL::squared_distance(c, q)) - r; +} + +template +typename K::FT +box(const typename K::Point_3& b, + const typename K::Point_3& t, + const typename K::Point_3& q) +{ + typename K::Point_3 c = CGAL::midpoint(b, t); + typename K::Iso_cuboid_3 ic(b, t); + bool inside = ic.has_on_bounded_side(q); + typename K::FT d = 0; + if(inside) + { + d = (std::min)({CGAL::abs(q.x() - b.x()), CGAL::abs(q.x() - t.x()), + CGAL::abs(q.y() - b.y()), CGAL::abs(q.y() - t.y()), + CGAL::abs(q.z() - b.z()), CGAL::abs(q.z() - t.z())}); + } + else + { + for(int i=0; i<3; ++i) + d += (CGAL::abs(q[i] - c[i]) > (c[i] - b[i]) ? CGAL::square(q[i] - c[i]) : 0); + d = CGAL::approximate_sqrt(d); + } + + return inside ? - d : d; +} + +// template +// typename K::FT +// disk(const typename K::Point_3& c, +// const typename K::Vector_3& n, +// const typename K::FT r, +// const typename K::Point_3& q) +// { +// typename K::Plane_3 pl(c, n); +// typename K::Point_3 pq = pl.projection(q); + +// typename K::FT sq_dpl = CGAL::squared_distance(q, pl); + +// if(CGAL::squared_distance(pq, c) < CGAL::square(r)) +// return CGAL::approximate_sqrt(sq_dpl); +// else +// return CGAL::approximate_sqrt(CGAL::square(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - r) + sq_dpl); +// } + + + +// p is the center of the base disk +// q is the center of the top disk +template +typename K::FT +infinite_cylinder(const typename K::Point_3& b, + const typename K::Vector_3& n, + const typename K::FT r, + const typename K::Point_3& q) +{ + typename K::Plane_3 pl(b, n); + typename K::Point_3 pq = pl.projection(q); + return CGAL::approximate_sqrt(CGAL::squared_distance(pq, b)) - r; +} + + +// c is the center of the torus +// n is the normal of the plane containing all centers of the tube +// r is the small radius +// R is the large radius +template +typename K::FT +torus(const typename K::Point_3& c, + const typename K::Vector_3& n, + const typename K::FT r, + const typename K::FT R, + const typename K::Point_3& q) +{ + typename K::Vector_3 w (c, q); + typename K::Plane_3 pl(c, n); + typename K::Point_3 pq = pl.projection(q); + typename K::FT d = CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R; + typename K::FT h = CGAL::abs(CGAL::scalar_product(w, n)); + + return CGAL::approximate_sqrt(CGAL::square(d) + CGAL::square(h)) - r; +} + +template +typename K::FT +torus_ridge(const typename K::Point_3& c, + const typename K::Vector_3& n, + const typename K::FT r, + const typename K::FT R, + const typename K::Point_3& q) +{ + typename K::Vector_3 w (c, q); + typename K::Plane_3 pl(c, n); + typename K::Point_3 pq = pl.projection(q); + typename K::FT d = CGAL::abs(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R) - r; + return d + CGAL::squared_distance(q, pl); +} + +template +typename K::FT +inverted_torus(const typename K::Point_3& c, + const typename K::Vector_3& n, + const typename K::FT r, + const typename K::FT R, + const typename K::Point_3& q) +{ + typename K::Vector_3 w (c, q); + typename K::Plane_3 pl(c, n); + typename K::Point_3 pq = pl.projection(q); + typename K::FT d = CGAL::abs(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - R) - r; + return d - CGAL::squared_distance(q, pl); +} + +///////////////////////////////////////////////////////////////// + +template +typename K::FT +shape_union(const S1& s1, const S2& s2, const typename K::Point_3& q) +{ + return std::min(s1(q), s2(q)); +} + +template +typename K::FT +shape_difference(const S1& s1, const S2& s2, const typename K::Point_3& q) +{ + return std::max(s1(q), -s2(q)); +} + +template +typename K::FT +shape_intersection(const S1& s1, const S2& s2, const typename K::Point_3& q) +{ + return std::max(s1(q), s2(q)); +} + +template +typename K::FT +shape_symmetric_difference(const S1& s1, const S2& s2, const typename K::Point_3& q) +{ + return std::max(-std::min(s1(q), s2(q)), std::max(s1(q), s2(q))); +} + + +} // namespace Shapes +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERNAL_IMPLICIT_SHAPES_HELPER_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h index 7090315ddfb..1ffc7da7fcf 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h @@ -76,7 +76,12 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi typename GeomTraits::Compute_z_3 z_coord = gt.compute_z_3_object(); typename GeomTraits::Construct_point_3 point = gt.construct_point_3_object(); - FT mu = FT(0.0); + FT mu = FT(0); + + // @todo, technically we should be using the edge intersection oracle here, but there is a nuance + // between MC and DC on the handling of edges that have val0 = val1 = isovalue: in MC we assume + // the isosurface is in the middle, in DC we assume the isosurface is not intersecting the edge. + // In the oracle, we follow DC right now. Could put a Boolean parameter, but it's ugly. // don't divide by 0 if(abs(d1 - d0) < 0.000001) // @fixme hardcoded bound @@ -87,9 +92,9 @@ typename GeomTraits::Point_3 vertex_interpolation(const typename GeomTraits::Poi CGAL_assertion(mu >= FT(0.0) || mu <= FT(1.0)); // linear interpolation - return point(x_coord(p1) * mu + x_coord(p0) * (FT(1.0) - mu), - y_coord(p1) * mu + y_coord(p0) * (FT(1.0) - mu), - z_coord(p1) * mu + z_coord(p0) * (FT(1.0) - mu)); + return point(x_coord(p1) * mu + x_coord(p0) * (FT(1) - mu), + y_coord(p1) * mu + y_coord(p0) * (FT(1) - mu), + z_coord(p1) * mu + z_coord(p0) * (FT(1) - mu)); } // retrieves the corner vertices and their values of a cell and return the lookup index @@ -128,7 +133,7 @@ template -void mc_construct_vertices(const typename Domain::Cell_descriptor cell, +void MC_construct_vertices(const typename Domain::Cell_descriptor cell, const std::size_t i_case, const Corners& corners, const Values& values, @@ -145,9 +150,9 @@ void mc_construct_vertices(const typename Domain::Cell_descriptor cell, std::size_t flag = 1; std::size_t e_id = 0; - for(const Edge_descriptor& edge : cell_edges) + for(const Edge_descriptor& e : cell_edges) { - (void)edge; // @todo + CGAL_USE(e); if(flag & Cube_table::intersected_edges[i_case]) { @@ -187,11 +192,11 @@ void mc_construct_triangles(const int i_case, const int eg2 = Cube_table::triangle_cases[t_index + 2]; // insert new triangle in list - #ifdef CGAL_LINKED_WITH_TBB +#ifdef CGAL_LINKED_WITH_TBB auto& tris = triangles.local(); - #else +#else auto& tris = triangles; - #endif +#endif tris.push_back({vertices[eg0], vertices[eg1], vertices[eg2]}); } } @@ -203,12 +208,12 @@ void triangles_to_polygon_soup(const TriangleRange& triangles, PointRange& points, PolygonRange& polygons) { - #ifdef CGAL_LINKED_WITH_TBB +#ifdef CGAL_LINKED_WITH_TBB for(const auto& triangle_list : triangles) { - #else +#else const auto& triangle_list = triangles; - #endif +#endif for(const auto& triangle : triangle_list) { @@ -222,9 +227,9 @@ void triangles_to_polygon_soup(const TriangleRange& triangles, polygons.push_back({id + 2, id + 1, id + 0}); } - #ifdef CGAL_LINKED_WITH_TBB +#ifdef CGAL_LINKED_WITH_TBB } - #endif +#endif } // Marching Cubes implemented as a functor that runs on every cell of the grid @@ -270,7 +275,6 @@ public: // computes one cell void operator()(const Cell_descriptor& cell) { - // @todo: maybe better checks if the domain can be processed? CGAL_precondition(m_domain.cell_vertices(cell).size() == 8); CGAL_precondition(m_domain.cell_edges(cell).size() == 12); @@ -284,7 +288,7 @@ public: return; std::array vertices; - mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices); + MC_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices); mc_construct_triangles(i_case, vertices, m_triangles); } diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits.h new file mode 100644 index 00000000000..46229fde32f --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits.h @@ -0,0 +1,26 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Mael Rouxel-Labbé + +#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H +#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H + +#include + +namespace CGAL { +namespace Isosurfacing { + +template +struct partition_traits; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Cartesian_grid.h similarity index 60% rename from Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h rename to Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Cartesian_grid.h index 0c991c4a4c0..455f2063c7d 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Cartesian_grid.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,16 +8,17 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Julian Stahl +// Mael Rouxel-Labbé -#ifndef CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H -#define CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H +#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H +#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H #include +#include #include #include -#include #include #ifdef CGAL_LINKED_WITH_TBB @@ -29,13 +30,18 @@ namespace CGAL { namespace Isosurfacing { -namespace internal { -// The topology of a Cartesian grid. -// All elements are created with the help of the cube tables. -class Grid_topology_3 +template +class Cartesian_grid_3; + +template +struct partition_traits; + +template +struct partition_traits > { -public: + using Self = Cartesian_grid_3; + // identifies a vertex by its (i, j, k) indices using Vertex_descriptor = std::array; @@ -54,20 +60,15 @@ public: using Cell_vertices = std::array; using Cell_edges = std::array; -public: - // creates the topology of a grid with size_i, size_j, and size_k vertices in the respective dimensions. - Grid_topology_3(const std::size_t size_i, - const std::size_t size_j, - const std::size_t size_k) - : size_i{size_i}, - size_j{size_j}, - size_k{size_k} + static decltype(auto) /*Point_3*/ point(const Vertex_descriptor& v, + const Self& g) { - CGAL_precondition(size_i > 0 && size_j > 0 && size_k > 0); + return g.point(v[0], v[1], v[2]); } // gets a container with the two vertices incident to edge e - Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e) const + static Vertices_incident_to_edge incident_vertices(const Edge_descriptor& e, + const Self&) { Vertices_incident_to_edge ev; ev[0] = { e[0], e[1], e[2] }; // start vertex @@ -77,7 +78,8 @@ public: } // gets a container with all cells incident to edge e - Cells_incident_to_edge incident_cells(const Edge_descriptor& e) const + static Cells_incident_to_edge incident_cells(const Edge_descriptor& e, + const Self&) { // lookup the neighbor cells relative to the edge const int local = internal::Cube_table::edge_store_index[e[3]]; @@ -96,7 +98,8 @@ public: } // gets a container with all vertices of cell c - Cell_vertices cell_vertices(const Cell_descriptor& c) const + static Cell_vertices cell_vertices(const Cell_descriptor& c, + const Self&) { Cell_vertices cv; for(std::size_t i=0; i - void for_each_vertex(Functor& f, Sequential_tag) const + static void for_each_vertex(Functor& f, + const Self& g, + const CGAL::Sequential_tag) { - for(std::size_t i=0; i - void for_each_edge(Functor& f, Sequential_tag) const + static void for_each_edge(Functor& f, + const Self& g, + const CGAL::Sequential_tag) { - for(std::size_t i=0; i - void for_each_cell(Functor& f, Sequential_tag) const + static void for_each_cell(Functor& f, + const Self& g, + const CGAL::Sequential_tag) { - for(std::size_t i=0; i - void for_each_vertex(Functor& f, Parallel_tag) const + static void for_each_vertex(Functor& f, + const Self& g, + const CGAL::Parallel_tag) { - const std::size_t sj = size_j; - const std::size_t sk = size_k; + const std::size_t sj = g.ydim(); + const std::size_t sk = g.zdim(); // for now only parallelize outer loop auto iterator = [&f, sj, sk](const tbb::blocked_range& r) { - for(std::size_t i = r.begin(); i != r.end(); ++i) + for(std::size_t i=r.begin(); i!=r.end(); ++i) for(std::size_t j=0; j(0, size_i), iterator); + tbb::parallel_for(tbb::blocked_range(0, g.xdim()), iterator); } // iterates in parallel over all edges e calling f(e) on every one template - void for_each_edge(Functor& f, Parallel_tag) const + static void for_each_edge(Functor& f, + const Self& g, + const CGAL::Parallel_tag) { - const std::size_t sj = size_j; - const std::size_t sk = size_k; + const std::size_t sj = g.ydim(); + const std::size_t sk = g.zdim(); // for now only parallelize outer loop auto iterator = [&f, sj, sk](const tbb::blocked_range& r) { - for(std::size_t i = r.begin(); i != r.end(); ++i) { - for(std::size_t j=0; j(0, size_i - 1), iterator); + tbb::parallel_for(tbb::blocked_range(0, g.xdim() - 1), iterator); } // iterates in parallel over all cells c calling f(c) on every one template - void for_each_cell(Functor& f, Parallel_tag) const + static void for_each_cell(Functor& f, + const Self& g, + const CGAL::Parallel_tag) { - const std::size_t sj = size_j; - const std::size_t sk = size_k; + const std::size_t sj = g.ydim(); + const std::size_t sk = g.zdim(); // for now only parallelize outer loop - auto iterator = [&f, sj, sk](const tbb::blocked_range3d& r) { + auto iterator = [&f, sj, sk](const tbb::blocked_range3d& r) + { const std::size_t i_begin = r.pages().begin(); const std::size_t i_end = r.pages().end(); const std::size_t j_begin = r.rows().begin(); @@ -232,19 +250,13 @@ public: f({i, j, k}); }; - tbb::blocked_range3d range(0, size_i - 1, 0, size_j - 1, 0, size_k - 1); + tbb::blocked_range3d range(0, g.xdim() - 1, 0, g.ydim() - 1, 0, g.zdim() - 1); tbb::parallel_for(range, iterator); } -#endif // CGAL_LINKED_WITH_TBB - -private: - std::size_t size_i; - std::size_t size_j; - std::size_t size_k; + #endif // CGAL_LINKED_WITH_TBB }; -} // namespace internal } // namespace Isosurfacing } // namespace CGAL -#endif // CGAL_ISOSURFACING_3_INTERNAL_GRID_TOPOLOGY_3_H +#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_CARTESIAN_GRID_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h new file mode 100644 index 00000000000..cdd69e353c8 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h @@ -0,0 +1,27 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Julian Stahl + +#ifndef CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H +#define CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H + +#include + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +// @todo + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERNAL_PARTITION_TRAITS_OCTREE_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/tables.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/tables.h index 6a39a273c3c..9b0771ffcf7 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/tables.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/tables.h @@ -46,7 +46,6 @@ namespace CGAL { namespace Isosurfacing { namespace internal { - namespace Cube_table { /* * Naming convention from "A parallel dual marching cubes approach @@ -681,7 +680,6 @@ constexpr int t_ambig[256] = }; } // namespace Cube_table - } // namespace internal } // namespace Isosurfacing } // namespace CGAL diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h index 9c6c4c44e2a..9031639b0b8 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h @@ -46,8 +46,10 @@ #include #include -#include -#include +#ifdef CGAL_LINKED_WITH_TBB +# include +# include +#endif #include #include @@ -94,24 +96,32 @@ private: } }; - using Edge_point_map = tbb::concurrent_hash_map; - private: const Domain& m_domain; FT m_isovalue; std::atomic m_point_counter; + +#ifdef CGAL_LINKED_WITH_TBB tbb::concurrent_vector m_points; - - Edge_point_map m_edges; - tbb::concurrent_vector > m_triangles; + using Edge_point_map = tbb::concurrent_hash_map; + Edge_point_map m_edges; +#else + std::vector m_points; + std::vector > m_triangles; + + // std::unordered_map m_edges; // @tmp hash map + std::map m_edges; // @tmp hash map +#endif + public: TMC_functor(const Domain& domain, const FT isovalue) : m_domain(domain), - m_isovalue(isovalue) + m_isovalue(isovalue), + m_point_counter(0) { } void operator()(const Cell_descriptor& cell) @@ -124,7 +134,7 @@ public: const int tcm = Cube_table::t_ambig[i_case]; if(tcm == 105) { - if (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; @@ -132,12 +142,10 @@ public: constexpr int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1 if(i_case == 0 || i_case == all_bits_set) - { return; - } std::array vertices; - mc_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices); + MC_construct_vertices(cell, i_case, corners, values, m_isovalue, m_domain, vertices); // @todo improve triangle generation @@ -195,28 +203,47 @@ private: bool find_point(const Edge_index& e, Point_index& i) { +#ifdef CGAL_LINKED_WITH_TBB typename Edge_point_map::const_accessor acc; if (m_edges.find(acc, e)) { i = acc->second; return true; } +#else + auto it = m_edges.find(e); + if (it != m_edges.end()) + { + i = it->second; + return true; + } +#endif return false; } Point_index add_point(const Point_3& p, const Edge_index& e) { + const Point_index i = m_point_counter++; + std::cout << "i =" << i << std::endl; + +#ifdef CGAL_LINKED_WITH_TBB typename 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; +#else + auto res = m_edges.insert({e, i}); + if (!res.second) + return res.first->second; + m_points.resize(i + 1); +#endif + + m_points[i] = p; return i; } @@ -224,7 +251,11 @@ private: { const Point_index i = m_point_counter++; +#ifdef CGAL_LINKED_WITH_TBB m_points.grow_to_at_least(i + 1); +#else + m_points.resize(i + 1); +#endif m_points[i] = p; return i; @@ -248,8 +279,6 @@ 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(); - using uint = unsigned int; - // 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}; auto get_edge_vertex = [](const int e, unsigned int& v0, unsigned int& v1, const unsigned char l_edges_[12]) @@ -272,7 +301,7 @@ private: if(flag & Cube_table::intersected_edges[i_case]) { // generate vertex here, do not care at this point if vertex already exists - uint v0, v1; + unsigned int v0, v1; get_edge_vertex(eg, v0, v1, l_edges_); FT l = (i0 - values[v0]) / (values[v1] - values[v0]); @@ -358,15 +387,15 @@ private: for(int f=0; f<6; ++f) { // classify face - unsigned int f_case{0}; - uint v0 = get_face_v(f, 0); - uint v1 = get_face_v(f, 1); - uint v2 = get_face_v(f, 2); - uint v3 = get_face_v(f, 3); - uint e0 = get_face_e(f, 0); - uint e1 = get_face_e(f, 1); - uint e2 = get_face_e(f, 2); - uint e3 = get_face_e(f, 3); + unsigned int f_case = 0; + unsigned int v0 = get_face_v(f, 0); + unsigned int v1 = get_face_v(f, 1); + unsigned int v2 = get_face_v(f, 2); + unsigned int v3 = get_face_v(f, 3); + unsigned int e0 = get_face_e(f, 0); + unsigned int e1 = get_face_e(f, 1); + unsigned int e2 = get_face_e(f, 2); + unsigned int e3 = get_face_e(f, 3); FT f0 = values[v0]; FT f1 = values[v1]; FT f2 = values[v2]; @@ -584,9 +613,9 @@ private: // set corresponging edge auto set_c = [](const int cnt, const int pos, const int val, unsigned long long& c_) { - const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; - const uint c_sz = c_ & mask[cnt]; - const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); + const unsigned int mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; + const unsigned int c_sz = c_ & mask[cnt]; + const unsigned int e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); c_ &= ~(((unsigned long long)0xF) << e); c_ |= (((unsigned long long)val) << e); }; @@ -594,22 +623,22 @@ private: // read edge from contour auto get_c = [](const int cnt, const int pos, unsigned long long c_) -> int { - const uint mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; - const uint c_sz = (uint)(c_ & mask[cnt]); - const uint e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); + const unsigned int mask[4] = {0x0, 0xF, 0xFF, 0xFFF}; + const unsigned int c_sz = (unsigned int)(c_ & mask[cnt]); + const unsigned int e = 16 + 4 * ((c_sz & 0xF) + ((c_sz & 0xF0) >> 4) + ((c_sz & 0xF00) >> 8) + pos); return int((c_ >> e) & 0xF); }; // connect oriented contours - uint cnt_{0}; - for(uint e=0; e<12; ++e) + unsigned int cnt_ = 0; + for(unsigned int e=0; e<12; ++e) { if(is_segm_set(e, segm_)) { - uint eTo = get_segm(e, 0, segm_); - uint eIn = get_segm(e, 1, segm_); - uint eStart = e; - uint pos = 0; + unsigned int eTo = get_segm(e, 0, segm_); + unsigned int eIn = get_segm(e, 1, segm_); + unsigned int eStart = e; + unsigned int pos = 0; set_c(cnt_, pos, eStart, c_); while(eTo != eStart) @@ -637,7 +666,7 @@ private: FT ui[2]{}; FT vi[2]{}; FT wi[2]{}; - unsigned char q_sol{0}; + unsigned char q_sol = 0; const FT a = (values[0] - values[1]) * (-values[6] + values[7] + values[4] - values[5]) - (values[4] - values[5]) * (-values[2] + values[3] + values[0] - values[1]); const FT b = (i0 - values[0]) * (-values[6] + values[7] + values[4] - values[5]) + @@ -747,16 +776,16 @@ private: auto numberOfSetBits = [](const unsigned char n) { // C or C++: use uint32_t - uint b = uint{n}; + unsigned int b = (unsigned int)(n); b = b - ((b >> 1) & 0x55555555); b = (b & 0x33333333) + ((b >> 2) & 0x33333333); return (((b + (b >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24; }; // compute the number of solutions to the quadratic equation for a given face - auto nrQSolFace = [](const uint f, const unsigned char n) + auto nrQSolFace = [](const unsigned int f, const unsigned char n) { - uint nr{0}; + unsigned int nr = 0; switch (f) { case 0: @@ -837,9 +866,9 @@ private: { FT umin(2); FT umax(-2); - const uint e0 = get_c(t, 0, c_); - const uint e1 = get_c(t, 1, c_); - const uint e2 = get_c(t, 2, c_); + const unsigned int e0 = get_c(t, 0, c_); + const unsigned int e1 = get_c(t, 1, c_); + const unsigned int e2 = get_c(t, 2, c_); const FT u_e0 = e_vert(e0, 0); const FT u_e1 = e_vert(e1, 0); const FT u_e2 = e_vert(e2, 0); @@ -894,9 +923,9 @@ private: const int cnt_sz = int(get_cnt_size(i, c_)); for(int r=0; r::max(); - uint ci = get_c(i, r, c_); + unsigned int ci = get_c(i, r, c_); const FT u_edge = e_vert(ci, 0); const FT v_edge = e_vert(ci, 1); const FT w_edge = e_vert(ci, 2); @@ -934,10 +963,11 @@ private: for(int r=0; r + +#include +#include + +namespace CGAL { +namespace Isosurfacing { + +/*! + * \ingroup IS_Fields_helpers_grp + * + * \cgalModels{InterpolationScheme_3} + * + * The class `Trilinear_interpolation` is the standard interpolation scheme to extrapolate + * data defined only at vertices of a Cartesian grid. + * + * \tparam Grid must be `CGAL::Isosurfacing::Cartesian_grid_3`, with `GeomTraits` + * a model of `IsosurfacingTraits_3` + */ +template +class Trilinear_interpolation +{ +public: + using Geom_traits = typename Grid::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; + using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3; + +public: + /*! + * \brief interpolates the values at a given point using trilinear interpolation. + * + * \param p the point at which to interpolate the values + * \param g the grid + * \param values the continuous field of scalar values, defined over the bounding box of `g` + * + * \return the interpolated value at point `p` + */ + FT interpolate_values(const Point_3& p, + const Grid& g, + const std::vector& values) const +{ + typename Geom_traits::Compute_x_3 x_coord = g.geom_traits().compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = g.geom_traits().compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = g.geom_traits().compute_z_3_object(); + typename Geom_traits::Construct_vertex_3 vertex = g.geom_traits().construct_vertex_3_object(); + + // trilinear interpolation of stored values + const Iso_cuboid_3& bbox = g.bbox(); + const std::array& spacing = g.spacing(); + + // calculate min index including border case + const Point_3& min_p = vertex(bbox, 0); + std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0]; + std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1]; + std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2]; + + // @todo check this + if(i == g.xdim() - 1) + --i; + if(j == g.ydim() - 1) + --j; + if(k == g.zdim() - 1) + --k; + + // calculate coordinates of min index + const FT min_x = i * spacing[0] + x_coord(min_p); + const FT min_y = j * spacing[1] + y_coord(min_p); + const FT min_z = k * spacing[2] + z_coord(min_p); + + // interpolation factors between 0 and 1 + const FT f_i = (x_coord(p) - min_x) / spacing[0]; + const FT f_j = (y_coord(p) - min_y) / spacing[1]; + const FT f_k = (z_coord(p) - min_z) / spacing[2]; + + // read the value at all 8 corner points + const FT g000 = values[g.linear_index(i + 0, j + 0, k + 0)]; + const FT g001 = values[g.linear_index(i + 0, j + 0, k + 1)]; + const FT g010 = values[g.linear_index(i + 0, j + 1, k + 0)]; + const FT g011 = values[g.linear_index(i + 0, j + 1, k + 1)]; + const FT g100 = values[g.linear_index(i + 1, j + 0, k + 0)]; + const FT g101 = values[g.linear_index(i + 1, j + 0, k + 1)]; + const FT g110 = values[g.linear_index(i + 1, j + 1, k + 0)]; + const FT g111 = values[g.linear_index(i + 1, j + 1, k + 1)]; + + // interpolate along all axes by weighting the corner points + const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k); + const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k; + const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k); + const FT lambda011 = (FT(1) - f_i) * f_j * f_k; + const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k); + const FT lambda101 = f_i * (FT(1) - f_j) * f_k; + const FT lambda110 = f_i * f_j * (FT(1) - f_k); + const FT lambda111 = f_i * f_j * f_k; + + // add weighted corners + return g000 * lambda000 + g001 * lambda001 + + g010 * lambda010 + g011 * lambda011 + + g100 * lambda100 + g101 * lambda101 + + g110 * lambda110 + g111 * lambda111; + } + + /*! + * \brief interpolates the gradients at a given point using trilinear interpolation. + * + * \param p the point at which to interpolate the gradients + * \param g the grid + * \param gradients the continuous field of vector values, defined over the bounding box of `g` + * + * \return the interpolated value at point `p` + */ + Vector_3 interpolate_gradients(const Point_3& p, + const Grid& g, + const std::vector& gradients) const + { + typename Geom_traits::Compute_x_3 x_coord = g.geom_traits().compute_x_3_object(); + typename Geom_traits::Compute_y_3 y_coord = g.geom_traits().compute_y_3_object(); + typename Geom_traits::Compute_z_3 z_coord = g.geom_traits().compute_z_3_object(); + typename Geom_traits::Construct_vector_3 vector = g.geom_traits().construct_vector_3_object(); + typename Geom_traits::Construct_vertex_3 vertex = g.geom_traits().construct_vertex_3_object(); + + // trilinear interpolation of stored gradients + const Iso_cuboid_3& bbox = g.bbox(); + const std::array& spacing = g.spacing(); + + // calculate min index including border case + const Point_3& min_p = vertex(bbox, 0); + std::size_t i = (x_coord(p) - x_coord(min_p)) / spacing[0]; + std::size_t j = (y_coord(p) - y_coord(min_p)) / spacing[1]; + std::size_t k = (z_coord(p) - z_coord(min_p)) / spacing[2]; + + if(i == g.xdim() - 1) + --i; + if(j == g.ydim() - 1) + --j; + if(k == g.zdim() - 1) + --k; + + // calculate coordinates of min index + const FT min_x = i * spacing[0] + x_coord(min_p); + const FT min_y = j * spacing[1] + y_coord(min_p); + const FT min_z = k * spacing[2] + z_coord(min_p); + + // interpolation factors between 0 and 1 + const FT f_i = (x_coord(p) - min_x) / spacing[0]; + const FT f_j = (y_coord(p) - min_y) / spacing[1]; + const FT f_k = (z_coord(p) - min_z) / spacing[2]; + + // read the value at all 8 corner points + const Vector_3& g000 = gradients[g.linear_index(i + 0, j + 0, k + 0)]; + const Vector_3& g001 = gradients[g.linear_index(i + 0, j + 0, k + 1)]; + const Vector_3& g010 = gradients[g.linear_index(i + 0, j + 1, k + 0)]; + const Vector_3& g011 = gradients[g.linear_index(i + 0, j + 1, k + 1)]; + const Vector_3& g100 = gradients[g.linear_index(i + 1, j + 0, k + 0)]; + const Vector_3& g101 = gradients[g.linear_index(i + 1, j + 0, k + 1)]; + const Vector_3& g110 = gradients[g.linear_index(i + 1, j + 1, k + 0)]; + const Vector_3& g111 = gradients[g.linear_index(i + 1, j + 1, k + 1)]; + + // interpolate along all axes by weighting the corner points + const FT lambda000 = (FT(1) - f_i) * (FT(1) - f_j) * (FT(1) - f_k); + const FT lambda001 = (FT(1) - f_i) * (FT(1) - f_j) * f_k; + const FT lambda010 = (FT(1) - f_i) * f_j * (FT(1) - f_k); + const FT lambda011 = (FT(1) - f_i) * f_j * f_k; + const FT lambda100 = f_i * (FT(1) - f_j) * (FT(1) - f_k); + const FT lambda101 = f_i * (FT(1) - f_j) * f_k; + const FT lambda110 = f_i * f_j * (FT(1) - f_k); + const FT lambda111 = f_i * f_j * f_k; + + // add weighted corners + return vector(x_coord(g000) * lambda000 + x_coord(g001) * lambda001 + + x_coord(g010) * lambda010 + x_coord(g011) * lambda011 + + x_coord(g100) * lambda100 + x_coord(g101) * lambda101 + + x_coord(g110) * lambda110 + x_coord(g111) * lambda111, + y_coord(g000) * lambda000 + y_coord(g001) * lambda001 + + y_coord(g010) * lambda010 + y_coord(g011) * lambda011 + + y_coord(g100) * lambda100 + y_coord(g101) * lambda101 + + y_coord(g110) * lambda110 + y_coord(g111) * lambda111, + z_coord(g000) * lambda000 + z_coord(g001) * lambda001 + + z_coord(g010) * lambda010 + z_coord(g011) * lambda011 + + z_coord(g100) * lambda100 + z_coord(g101) * lambda101 + + z_coord(g110) * lambda110 + z_coord(g111) * lambda111); + } +}; + +// This can be used for example when we have implicit functions for data (values & gradients), +// but use an interpolated values field as to store data. +template +class Function_evaluation +{ + using Geom_traits = typename Grid::Geom_traits; + using FT = typename Geom_traits::FT; + using Point_3 = typename Geom_traits::Point_3; + using Vector_3 = typename Geom_traits::Vector_3; + + std::function m_value_fn; + std::function m_gradient_fn; + +public: + template + Function_evaluation(const ValueFunction& value_fn, + const GradientFunction& gradient_fn = [](const Point_3&) -> Vector_3 { return CGAL::NULL_VECTOR; }) + : m_value_fn{value_fn}, + m_gradient_fn{gradient_fn} + { } + +public: + FT interpolate_values(const Point_3& p, const Grid&, const std::vector&) const + { + return m_value_fn(p); + } + + Vector_3 interpolate_gradients(const Point_3& p, const Grid&, const std::vector&) const + { + return m_gradient_fn(p); + } +}; + +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERNAL_INTERPOLATION_SCHEMES_3_H diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h index cd0a22520c2..c84808be4e7 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2022-2023 INRIA Sophia-Antipolis (France). +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,6 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Julian Stahl +// Mael Rouxel-Labbé #ifndef CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H #define CGAL_ISOSURFACING_3_MARCHING_CUBES_3_H @@ -38,20 +39,21 @@ namespace Isosurfacing { * and `BackInsertionSequence` whose value type is `std::size_t`. * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * - * \param domain the domain providing input data and its topology - * \param isovalue value of the isosurface - * \param points points of the triangles in the created indexed face set + * \param domain the domain providing the spacial partition and the data + * \param isovalue the value defining the isosurface + * \param points the 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 np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below * * \cgalNamedParamsBegin * \cgalParamNBegin{use_topologically_correct_marching_cubes} - * \cgalParamDescription{whether the topologically correct variant of Marching Cubes should be used} + * \cgalParamDescription{whether the topologically correct variant of Marching Cubes \cgalCite{cgal:g-ctcmi-16} should be used.} * \cgalParamType{Boolean} * \cgalParamDefault{`true`} * \cgalParamNEnd * \cgalNamedParamsEnd * + * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` */ template functor(domain, isovalue); domain.template for_each_cell(functor); functor.to_triangle_soup(points, triangles);