diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Orthtree_domain_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Orthtree_domain_3.h new file mode 100644 index 00000000000..7202a284018 --- /dev/null +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Orthtree_domain_3.h @@ -0,0 +1,182 @@ +// Copyright (c) 2022-2024 INRIA Sophia-Antipolis (France), GeometryFactory (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_INTERNAL_ORTHTREE_DOMAIN_3_H +#define CGAL_ISOSURFACING_3_INTERNAL_ORTHTREE_DOMAIN_3_H + +#include + +#include +#include +#include + +#include + +#include +#include + +namespace CGAL { +namespace Isosurfacing { +namespace internal { + +// 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 +// typedef Kernel_traits::Kernel::Triangle_3 Datum; +template +class Isosurfacing_domain_3::Kernel, std::vector >, ValueField, GradientField, IntersectionOracle> +{ +public: + using Edge_intersection_oracle = IntersectionOracle; + + using Partition = CGAL::Octree::Kernel, std::vector >; + + 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 PT = CGAL::Isosurfacing::partition_traits; + + using vertex_descriptor = typename PT::vertex_descriptor; + using edge_descriptor = typename PT::edge_descriptor; + using cell_descriptor = typename PT::cell_descriptor; + + using Edge_vertices = typename PT::Edge_vertices; + 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 = 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 Partition& m_partition; + const ValueField& m_values; + const GradientField& m_gradients; + const IntersectionOracle m_intersection_oracle; + mutable std::vector m_leaf_edges; // cache variable + mutable std::vector m_leaf_cells; // cache variable + +public: + 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_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. + + // returns the location of vertex `v` + decltype(auto) /*Point_3*/ point(const vertex_descriptor& v) const + { + return PT::point(v, m_partition); + } + + // returns the value of the function at vertex `v` + decltype(auto) /*FT*/ value(const vertex_descriptor& v) const + { + return m_values(v); + } + + // returns the value of the function at point `p` + decltype(auto) /*FT*/ value(const Point_3& p) const + { + return m_values(p); + } + + // returns the gradient at point `p` + decltype(auto) /*Vector_3*/ gradient(const Point_3& p) const + { + return m_gradients(p); + } + + // returns a container with the two vertices incident to the edge `e` + decltype(auto) /*Edge_vertices*/ incident_vertices(const edge_descriptor& e) const + { + return PT::incident_vertices(e, m_partition); + } + + // returns a container with all cells incident to the edge `e` + decltype(auto) /*Cells_incident_to_edge*/ incident_cells(const edge_descriptor& e) const + { + return PT::incident_cells(e, m_partition); + } + + // returns a container with all vertices of the cell `c` + decltype(auto) /*Cell_vertices*/ cell_vertices(const cell_descriptor& c) const + { + return PT::cell_vertices(c, m_partition); + } + + // returns a container with all edges of the cell `c` + decltype(auto) /*Cell_edges*/ cell_edges(const cell_descriptor& c) const + { + 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 + { + PT::template for_each_vertex(f, m_partition); + } + + // iterates over all edges `e`, calling `f(e)` on each of them + template + void for_each_edge(Functor& f) const + { + PT::template for_each_edge(f, m_leaf_edges, m_partition); + } + + // iterates over all cells `c`, calling `f(c)` on each of them + template + void for_each_cell(Functor& f) const + { + PT::template for_each_cell(f, m_leaf_cells, m_partition); + } + + // finds the intersection of the isosurface with the edge `e` (if any) + bool construct_intersection(const Point_3& p_0, const Point_3& p_1, + const FT val_0, const FT val_1, + const FT isovalue, + Point_3& p) const + { + return m_intersection_oracle(p_0, p_1, val_0, val_1, *this, isovalue, p); + } +}; + +} // namespace internal +} // namespace Isosurfacing +} // namespace CGAL + +#endif // CGAL_ISOSURFACING_3_INTERNAL_ORTHTREE_DOMAIN_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 index 2a4d1906a95..772790966e4 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/partition_traits_Octree.h @@ -263,14 +263,26 @@ public: const Orthtree& o, Sequential_tag) { - //for (const cell_descriptor& v : o.traverse(CGAL::Orthtrees::Leaves_traversal(o))) - // f(v); - std::set edge_set = get_leaf_edges(o); for(const edge_descriptor& e : o.leaf_edges()) f(e); } + template + static void for_each_edge(Functor& f, + std::vector& edges, + const Orthtree& o, + Sequential_tag) + { + if (edges.empty()) { + std::set edge_set = get_leaf_edges(o, edges); + std::copy(edge_set.begin(), edge_set.end(), std::back_inserter(edges)); + } + + for (const edge_descriptor& e : edges) + f(e); + } + template static void for_each_cell(Functor& f, const Orthtree& o, @@ -280,6 +292,20 @@ public: f(v); } + template + static void for_each_cell(Functor& f, + std::vector& cells, + const Orthtree& o, + CGAL::Sequential_tag) + { + if (cells.empty()) { + auto cell_range = o.traverse(CGAL::Orthtrees::Leaves_traversal(o)); + std::copy(cell_range.begin(), cell_range.end(), std::back_inserter(cells)); + } + for (const cell_descriptor& v : cells) + f(v); + } + #ifdef CGAL_LINKED_WITH_TBB template static void for_each_vertex(Functor& f, @@ -306,6 +332,25 @@ public: tbb::parallel_for_each(edges.begin(), edges.end(), f); } + template + static void for_each_edge(Functor& f, + std::vector& edges, + const Orthtree& o, + Parallel_tag) + { + if (edges.empty()) { + std::set edge_set = get_leaf_edges(o); + std::copy(edge_set.begin(), edge_set.end(), std::back_inserter(edges)); + } + + auto iterator = [&](const tbb::blocked_range& r) + { + for (std::size_t i = r.begin(); i != r.end(); ++i) + f(edges[i]); + }; + + tbb::parallel_for(tbb::blocked_range(0, edges.size()), iterator); + } template static void for_each_cell(Functor& f, @@ -321,13 +366,38 @@ public: const auto& cells = o.traverse(CGAL::Orthtrees::Leaves_traversal(o)); tbb::parallel_for_each(cells.begin(), cells.end(), func); } + template + static void for_each_cell(Functor& f, + std::vector& cells, + const Orthtree& o, + Parallel_tag) + { + if (cells.empty()) { + auto cell_range = o.traverse(CGAL::Orthtrees::Leaves_traversal(o)); + std::copy(cell_range.begin(), cell_range.end(), std::back_inserter(cells)); + } + + auto iterator = [&](const tbb::blocked_range& r) + { + for (std::size_t i = r.begin(); i != r.end(); ++i) { + const Uniform_coords& uc = uniform_coordinates(cells[i], o); + f(lex_index(uc[0], uc[1], uc[2], o.depth())); + } + }; + + tbb::parallel_for(tbb::blocked_range(0, cells.size()), iterator); + } #endif // CGAL_LINKED_WITH_TBB template static void for_each_vertex(Functor& f, const Orthtree& o) { return for_each_vertex(f, o, ConcurrencyTag{}); } template + static void for_each_edge(Functor& f, std::vector& edges, const Orthtree& o) { return for_each_edge(f, edges, o, ConcurrencyTag{}); } + template static void for_each_edge(Functor& f, const Orthtree& o) { return for_each_edge(f, o, ConcurrencyTag{}); } template + static void for_each_cell(Functor& f, std::vector& cells, const Orthtree& o) { return for_each_cell(f, cells, o, ConcurrencyTag{}); } + template static void for_each_cell(Functor& f, const Orthtree& o) { return for_each_cell(f, o, ConcurrencyTag{}); } private: diff --git a/Isosurfacing_3/test/Isosurfacing_3/dual_contouring_octree.cpp b/Isosurfacing_3/test/Isosurfacing_3/dual_contouring_octree.cpp index 41319a03203..5583d830259 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/dual_contouring_octree.cpp +++ b/Isosurfacing_3/test/Isosurfacing_3/dual_contouring_octree.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include