WIP adding Isosurfacing_domain_3 for Orthtree

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 <CGAL/license/Isosurfacing_3.h>
+
+#include <CGAL/Isosurfacing_3/internal/Cell_type.h>
+#include <CGAL/Isosurfacing_3/internal/partition_traits.h>
+#include <CGAL/Isosurfacing_3/edge_intersection_oracles_3.h>
+
+#include <CGAL/Octree.h>
+
+#include <CGAL/tags.h>
+#include <CGAL/Kernel_traits.h>
+
+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<Point>::Kernel::Triangle_3 Datum;
+template <typename ValueField,
+          typename GradientField,
+          typename IntersectionOracle>
+class Isosurfacing_domain_3<CGAL::Octree<typename Kernel_traits<typename ValueField::Point_3>::Kernel, std::vector<typename ValueField::Point_3> >, ValueField, GradientField, IntersectionOracle>
+{
+public:
+  using Edge_intersection_oracle = IntersectionOracle;
+
+  using Partition = CGAL::Octree<typename Kernel_traits<typename ValueField::Point_3>::Kernel, std::vector<typename ValueField::Point_3> >;
+
+  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<Partition>;
+
+  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<edge_descriptor> m_leaf_edges; // cache variable
+  mutable std::vector<cell_descriptor> 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 <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
+  void for_each_vertex(Functor& f) const
+  {
+    PT::template for_each_vertex<ConcurrencyTag>(f, m_partition);
+  }
+
+  // iterates over all edges `e`, calling `f(e)` on each of them
+  template <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
+  void for_each_edge(Functor& f) const
+  {
+    PT::template for_each_edge<ConcurrencyTag>(f, m_leaf_edges, m_partition);
+  }
+
+  // iterates over all cells `c`, calling `f(c)` on each of them
+  template <typename ConcurrencyTag = CGAL::Sequential_tag, typename Functor>
+  void for_each_cell(Functor& f) const
+  {
+    PT::template for_each_cell<ConcurrencyTag>(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

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<Orthtree>(o)))
-    //  f(v);
-
     std::set<edge_descriptor> edge_set = get_leaf_edges(o);
     for(const edge_descriptor& e : o.leaf_edges())
        f(e);
   }
 
+  template <typename Functor>
+  static void for_each_edge(Functor& f,
+                            std::vector<edge_descriptor>& edges,
+                            const Orthtree& o,
+                            Sequential_tag)
+  {
+    if (edges.empty()) {
+      std::set<edge_descriptor> 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 <typename Functor>
   static void for_each_cell(Functor& f,
                             const Orthtree& o,
@@ -280,6 +292,20 @@ public:
       f(v);
   }
 
+  template <typename Functor>
+  static void for_each_cell(Functor& f,
+                            std::vector<cell_descriptor>& cells,
+                            const Orthtree& o,
+                            CGAL::Sequential_tag)
+  {
+    if (cells.empty()) {
+      auto cell_range = o.traverse(CGAL::Orthtrees::Leaves_traversal<typename Orthtree::Orthtree>(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 <typename Functor>
   static void for_each_vertex(Functor& f,
@@ -306,6 +332,25 @@ public:
 
     tbb::parallel_for_each(edges.begin(), edges.end(), f);
   }
+  template <typename Functor>
+  static void for_each_edge(Functor& f,
+                            std::vector<edge_descriptor>& edges,
+                            const Orthtree& o,
+                            Parallel_tag)
+  {
+    if (edges.empty()) {
+      std::set<edge_descriptor> edge_set = get_leaf_edges(o);
+      std::copy(edge_set.begin(), edge_set.end(), std::back_inserter(edges));
+    }
+
+    auto iterator = [&](const tbb::blocked_range<std::size_t>& r)
+      {
+        for (std::size_t i = r.begin(); i != r.end(); ++i)
+          f(edges[i]);
+      };
+
+    tbb::parallel_for(tbb::blocked_range<std::size_t>(0, edges.size()), iterator);
+  }
 
   template <typename Functor>
   static void for_each_cell(Functor& f,
@@ -321,13 +366,38 @@ public:
     const auto& cells = o.traverse(CGAL::Orthtrees::Leaves_traversal<Orthtree>(o));
     tbb::parallel_for_each(cells.begin(), cells.end(), func);
   }
+  template <typename Functor>
+  static void for_each_cell(Functor& f,
+                            std::vector<cell_descriptor>& cells,
+                            const Orthtree& o,
+                            Parallel_tag)
+  {
+    if (cells.empty()) {
+      auto cell_range = o.traverse(CGAL::Orthtrees::Leaves_traversal<typename Orthtree::Orthtree>(o));
+      std::copy(cell_range.begin(), cell_range.end(), std::back_inserter(cells));
+    }
+
+    auto iterator = [&](const tbb::blocked_range<std::size_t>& 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<std::size_t>(0, cells.size()), iterator);
+  }
 #endif // CGAL_LINKED_WITH_TBB
 
   template <typename ConcurrencyTag, typename Functor>
   static void for_each_vertex(Functor& f, const Orthtree& o) { return for_each_vertex(f, o, ConcurrencyTag{}); }
   template <typename ConcurrencyTag, typename Functor>
+  static void for_each_edge(Functor& f, std::vector<edge_descriptor>& edges, const Orthtree& o) { return for_each_edge(f, edges, o, ConcurrencyTag{}); }
+  template <typename ConcurrencyTag, typename Functor>
   static void for_each_edge(Functor& f, const Orthtree& o) { return for_each_edge(f, o, ConcurrencyTag{}); }
   template <typename ConcurrencyTag, typename Functor>
+  static void for_each_cell(Functor& f, std::vector<cell_descriptor>& cells, const Orthtree& o) { return for_each_cell(f, cells, o, ConcurrencyTag{}); }
+  template <typename ConcurrencyTag, typename Functor>
   static void for_each_cell(Functor& f, const Orthtree& o) { return for_each_cell(f, o, ConcurrencyTag{}); }
 
 private:

Isosurfacing_3/test/Isosurfacing_3/dual_contouring_octree.cpp

@@ -6,6 +6,7 @@
 #include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
 #include <CGAL/Isosurfacing_3/Value_function_3.h>
 #include <CGAL/Isosurfacing_3/Gradient_function_3.h>
+#include <CGAL/Isosurfacing_3/internal/Orthtree_domain_3.h>
 #include <CGAL/Octree.h>
 
 #include <CGAL/IO/polygon_soup_io.h>