diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 09bde261d24b..8ec9ae609a36 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -38,6 +38,12 @@ Release date: October 2023 - Removed the class templates `Gray_image_mesh_domain_3`, `Implicit_mesh_domain_3`, and `Labeled_image_mesh_domain_3` which are deprecated since CGAL-4.13. +### [Quadtrees, Octrees, and Orthtrees](https://doc.cgal.org/6.0/Manual/packages.html#PkgOrthtree) +- **Breaking change**: + - Node splitting behavior and per-node data are now customizable via the Traits class. + - Nodes are now stored as a property map, with properties of each node accessed by index. + - Nearest neighbors functions only work for Orthtrees which provide the necessary functionality. + ### [Polygon Mesh Processing](https://doc.cgal.org/6.0/Manual/packages.html#PkgPolygonMeshProcessing) - Added the function `CGAL::Polygon_mesh_processing::interpolated_corrected_curvatures()` which can be used to compute diff --git a/Orthtree/benchmark/Orthtree/construction.cpp b/Orthtree/benchmark/Orthtree/construction.cpp index 36b8800f28c4..5d55f39c1ee8 100644 --- a/Orthtree/benchmark/Orthtree/construction.cpp +++ b/Orthtree/benchmark/Orthtree/construction.cpp @@ -12,16 +12,14 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; - -typedef CGAL::Search_traits_3 Kd_tree_traits; -typedef CGAL::Orthogonal_k_neighbor_search Kd_tree_search; -typedef Kd_tree_search::Tree Kdtree; +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; +using Kd_tree_traits = CGAL::Search_traits_3; +using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search; +using Kdtree = Kd_tree_search::Tree; using std::chrono::high_resolution_clock; using std::chrono::duration_cast; diff --git a/Orthtree/benchmark/Orthtree/nearest_neighbor.cpp b/Orthtree/benchmark/Orthtree/nearest_neighbor.cpp index 69ca01c89580..219d9f1b0f7a 100644 --- a/Orthtree/benchmark/Orthtree/nearest_neighbor.cpp +++ b/Orthtree/benchmark/Orthtree/nearest_neighbor.cpp @@ -12,21 +12,19 @@ #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; - -typedef CGAL::Search_traits_3 Kd_tree_traits; -typedef CGAL::Orthogonal_k_neighbor_search Kd_tree_search; -typedef Kd_tree_search::Tree Kdtree; - using std::chrono::high_resolution_clock; using std::chrono::duration_cast; using std::chrono::microseconds; +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; +using Kd_tree_traits = CGAL::Search_traits_3; +using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search; +using Kdtree = Kd_tree_search::Tree; + int main(int argc, char **argv) { int num_runs = 100; @@ -103,8 +101,8 @@ int main(int argc, char **argv) { // Time how long it takes to find neighbors using the octree auto octreeTime = bench( [&] { - std::vector nearest_neighbors; - octree.nearest_neighbors(search_point, k, std::back_inserter(nearest_neighbors)); + std::vector nearest_neighbors; + octree.nearest_k_neighbors(search_point, k, std::back_inserter(nearest_neighbors)); } ); diff --git a/Orthtree/benchmark/Orthtree/util.h b/Orthtree/benchmark/Orthtree/util.h index de65aa89e812..f619317fa5f5 100644 --- a/Orthtree/benchmark/Orthtree/util.h +++ b/Orthtree/benchmark/Orthtree/util.h @@ -10,8 +10,8 @@ template CGAL::Point_set_3 generate(size_t num_points = 1) { - typedef typename Kernel::Point_3 Point; - typedef CGAL::Point_set_3 Point_set; + using Point = typename Kernel::Point_3; + using Point_set = CGAL::Point_set_3; // Create an empty point set Point_set points; diff --git a/Orthtree/doc/Orthtree/Concepts/CollectionPartitioningOrthtreeTraits.h b/Orthtree/doc/Orthtree/Concepts/CollectionPartitioningOrthtreeTraits.h new file mode 100644 index 000000000000..97068f40a978 --- /dev/null +++ b/Orthtree/doc/Orthtree/Concepts/CollectionPartitioningOrthtreeTraits.h @@ -0,0 +1,101 @@ +/*! + \ingroup PkgOrthtreeConcepts + \cgalConcept + + Refinement of the `OrthtreeTraitsWithData` concept, adding requirements for the + traits class of a `CGAL::Orthtree` in order to supports nearest-neighbor searching. + + Nearest neighbor searches expect a tree where `Node_data` is a model of `ForwardRange`. + The leaf nodes of the tree represent an exclusive partition of the elements contained in the tree. + This means that no element should be contained by more than one node. + + \cgalRefines{OrthtreeTraitsWithData} + + \cgalHasModelsBegin + \cgalHasModels{CGAL::Orthtree_traits_point} + \cgalHasModelsEnd +*/ +class CollectionPartitioningOrthtreeTraits { +public: + + + /// \name Types + /// @{ + + /*! + * Sphere Type used for the shrinking-sphere approach for neighbor queries; needs to be copy assignable. + */ + using Sphere_d = unspecified_type; + + /*! + * \brief The data type contained by each node; must be a model of `ForwardRange` in addition to default constructible, copy constructible and copy assignable. + */ + using Node_data = unspecified_type; + + /*! + * \brief An element of the `Node_data` list-like type. + * + * Must be constructible from the value type of a `Node_data::iterator`. + * Typically the same as that type, but can also be an `std::reference_wrapper<>` if the type is not copyable. + */ + using Node_data_element = unspecified_type; + + /*! + * \brief Functor with an operator that calculates the squared distance of a `Node_data_element` from a point. + * + * Provides the operator: + * `FT operator()(const Node_data_element&, const Point_d&)` + */ + using Squared_distance_of_element = unspecified_type; + + /*! + * \brief Functor with an operator that constructs a `Sphere_d` from a provided center and squared radius. + * + * Provides the operator: + * `Sphere_d operator()(const Point_d&, const FT&)` + */ + using Construct_sphere_d = unspecified_type; + + /*! + * \brief Functor with an operator that provides the center of a `Sphere_d`. + * + * Provides the operator: + * `Point_d operator()(const Sphere_d&)` + */ + using Construct_center_d = unspecified_type; + + /*! + * \brief Functor with an operator that provides the squared radius of a `Sphere_d`. + * + * Provides the operator: + * `FT operator()(const Sphere_d&)` + */ + using Compute_squared_radius_d = unspecified_type; + + /// @} + + /// \name Operations + /// @{ + + /*! + * constructs an object of type `ConstructSphere_d`. + */ + Construct_sphere_d construct_sphere_d_object() const; + + /*! + * constructs an object of type `ConstructCenter_d`. + */ + Construct_center_d construct_center_d_object() const; + + /*! + * constructs an object of type `ComputeSquaredRadius_d`. + */ + Compute_squared_radius_d compute_squared_radius_d_object() const; + + /*! + * constructs an object of type `Squared_distance_of_element`. + */ + Squared_distance_of_element squared_distance_of_element_object() const; + + /// @} +}; diff --git a/Orthtree/doc/Orthtree/Concepts/OrthtreeTraits.h b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraits.h index fdb0781ba260..9f2e62d5db4a 100644 --- a/Orthtree/doc/Orthtree/Concepts/OrthtreeTraits.h +++ b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraits.h @@ -6,9 +6,9 @@ template parameter of the `CGAL::Orthtree` class. \cgalHasModelsBegin - \cgalHasModels{CGAL::Orthtree_traits_2} - \cgalHasModels{CGAL::Orthtree_traits_3} - \cgalHasModels{CGAL::Orthtree_traits_d} + \cgalHasModels{CGAL::Orthtree_traits_point} + \cgalHasModels{CGAL::Orthtree_traits_face_graph} + \cgalHasModels{CGAL::Orthtree_traits_base} \cgalHasModelsEnd */ class OrthtreeTraits @@ -17,31 +17,47 @@ class OrthtreeTraits /// \name Types /// @{ - - typedef unspecified_type Dimension; ///< Dimension type (see `CGAL::Dimension_tag`). - typedef unspecified_type Bbox_d; ///< Bounding box type. - typedef unspecified_type FT; ///< The number type of the %Cartesian coordinates of types `Point_d` - typedef unspecified_type Point_d; ///< Point type. - typedef unspecified_type Sphere_d; ///< The sphere type for neighbor queries. + using Node_index = unspecified_type; ///< An integer type for nodes + constexpr int dimension; ///< Dimension. + using FT = unspecified_type; ///< The number type of the %Cartesian coordinates of types `Point_d` + using Point_d = unspecified_type; ///< Point type. + using Bbox_d = unspecified_type; ///< Bounding box type. Must be constructible from a pair of `Point_d` objects. /*! A random access iterator type to enumerate the - %Cartesian coordinates of a point. + %Cartesian coordinates of a point of type `Point_d`. */ - typedef unspecified_type Cartesian_const_iterator_d; - typedef std::array Array; ///< Array used for easy point constructions. + using Cartesian_const_iterator_d = unspecified_type; - typedef unspecified_type Adjacency; ///< Specify the adjacency directions + /*! + * \brief Integral number type which can take on values indicating adjacency directions. + * + * Must be able to take on values ranging from 0 to the number of faces of the (hyper)rectangle type, equivalent to 2 * D. + */ + using Adjacency = unspecified_type; ///< Specify the adjacency directions /*! - Functor with an operator to construct a `Point_d` from an `Array` object. - */ - typedef unspecified_type Construct_point_d_from_array; + * \brief Functor with an operator to create the bounding box of the root node. + * + * Provides the operator: + * `Bbox_d operator()()` + * + * The bounding box must enclose all elements contained by the tree. + * It may be tight-fitting. The orthtree constructor produces a bounding box surrounding this region. + * For traits which assign no data to each node, this can be defined to return a fixed region. + */ + using Construct_root_node_bbox = unspecified_type; /*! - Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points). - */ - typedef unspecified_type Construct_bbox_d; + * \brief Functor with an operator to construct a `Point_d` from an initializer list of type `FT`. + * + * Provides the operator: + * `Point_d operator()(arg1, arg2,...)` + * + * For trees which use a different kernel for the bounding box type, + * the return type of this functor must match the kernel used by the bounding box type and not that of the contents. + */ + using Construct_point_d = unspecified_type; /// @} @@ -49,14 +65,14 @@ class OrthtreeTraits /// @{ /*! - Function used to construct an object of type `Construct_point_d_from_array`. - */ - Construct_point_d_from_array construct_point_d_from_array_object() const; + * constructs an object of type `Construct_root_node_bbox`. + */ + Construct_root_node_bbox construct_root_node_bbox_object() const; /*! - Function used to construct an object of type `Construct_bbox_d`. - */ - Construct_bbox_d construct_bbox_d_object() const; + * constructs an object of type `Construct_point_d`. + */ + Construct_point_d construct_point_d_object() const; /// @} }; diff --git a/Orthtree/doc/Orthtree/Concepts/OrthtreeTraitsWithData.h b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraitsWithData.h new file mode 100644 index 000000000000..f6533340079a --- /dev/null +++ b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraitsWithData.h @@ -0,0 +1,76 @@ +/*! + \ingroup PkgOrthtreeConcepts + \cgalConcept + + The concept `OrthtreeTraitsWithData` defines the requirements for the + template parameter of the `CGAL::Orthtree` class for a node type that stores data. + + \cgalRefines{OrthtreeTraits} + + \cgalHasModelsBegin + \cgalHasModels{CGAL::Orthtree_traits_point} + \cgalHasModels{CGAL::Orthtree_traits_face_graph} + \cgalHasModels{CGAL::Orthtree_traits_base} + \cgalHasModelsEnd +*/ +class OrthtreeTraitsWithData +{ +public: + + /// \name Types + /// @{ + /*! + * \brief The data type contained by each node. Must be default constructible, copy constructible and copy assignable. + */ + using Node_data = unspecified_type; + + /*! + * \brief Functor which initializes elements contained by the root node. + * + * Each node of a tree has an associated `Node_data` value. + * This functor initializes the `Node_data` of the root node. + * It takes no arguments, and returns an instance of `Node_data`. + * + * Provides the operator: + * `Node_data operator()()` + * + * Typically, the `Node_data` of the root node contains all the elements in the tree. + * For a tree in which each node contains a span (such as `std::span()`) this function would return the span containing all items. + * + */ + using Construct_root_node_contents = unspecified_type; + + /*! + * \brief Functor which fills the contents of the nodes children. + * + * Provides the operator: + * `void operator()(Node_index, Orthtree&, const Point_d&)` + * + * The functor is called during refinement of the `Orthtree` on a node after it has been split. The purpose of the functor is to + * distribute the `Node_data`, accessible via `tree.data()`, to the data of the nodes children, accessible via `tree.children()`. + * The first parameter is the `Node_index` of the node. The second parameter provides the instance of the `Orthtree` + * and the last parameter is the barycenter of the node which will be used as shared corner amongst the children of the node. + * + * For a tree in which each node contains a span, this may mean rearranging the contents of the original node + * and producing spans containing a subset of its contents for each of its children. + * For compatibility with locate, the center of the node is considered to be part of the upper half. + */ + using Distribute_node_contents = unspecified_type; + + /// @} + + /// \name Operations + /// @{ + + /*! + * constructs an object of type `Construct_root_node_contents`. + */ + Construct_root_node_contents construct_root_node_contents_object() const; + + /*! + * constructs an object of type `Distribute_node_contents`. + */ + Distribute_node_contents distribute_node_contents_object() const; + + /// @} +}; diff --git a/Orthtree/doc/Orthtree/Concepts/OrthtreeTraversal.h b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraversal.h index 5aeda690079e..64e84572cea5 100644 --- a/Orthtree/doc/Orthtree/Concepts/OrthtreeTraversal.h +++ b/Orthtree/doc/Orthtree/Concepts/OrthtreeTraversal.h @@ -1,13 +1,8 @@ - /*! \ingroup PkgOrthtreeConcepts \cgalConcept - \brief a traversal provides the functions needed to traverse the - nodes of an orthtree. - - A traversal is used to iterate on a tree with a user-selected order - (e.g. preorder, postorder). + \brief Requirements for defining a traversal strategy of an orthtree. \cgalHasModelsBegin \cgalHasModels{CGAL::Orthtrees::Preorder_traversal} @@ -20,15 +15,15 @@ class OrthtreeTraversal { public: - using Node = unspecified_type; ///< A specialization of [CGAL::Orthtree::Node](@ref CGAL::Orthtree::Node) + using Node_index = unspecified_type; ///< Index type of the orthtree to be traversed /*! - \brief returns the first node to iterate to, given the root of the Orthtree. + \brief returns the first node of the traversal. */ - Node first (Node root) const; + Node_index first_index() const; /*! - \brief returns the next node to iterate to, given the previous node. + \brief returns the next node to be traversed given the previous node `n`. */ - Node next(Node n) const; + Node_index next(Node_index n) const; }; diff --git a/Orthtree/doc/Orthtree/Doxyfile.in b/Orthtree/doc/Orthtree/Doxyfile.in index 980e12e8f6e5..2adfd640f365 100644 --- a/Orthtree/doc/Orthtree/Doxyfile.in +++ b/Orthtree/doc/Orthtree/Doxyfile.in @@ -5,3 +5,7 @@ PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Quadtrees, Octrees, and Orthtrees" EXTRACT_ALL = false HIDE_UNDOC_MEMBERS = true HIDE_UNDOC_CLASSES = true +WARN_IF_UNDOCUMENTED = false + +HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/orthtree.png \ + ${CGAL_PACKAGE_DOC_DIR}/fig/orthtree_mesh.png diff --git a/Orthtree/doc/Orthtree/Orthtree.txt b/Orthtree/doc/Orthtree/Orthtree.txt index bfe724038666..eb7b2559b2ac 100644 --- a/Orthtree/doc/Orthtree/Orthtree.txt +++ b/Orthtree/doc/Orthtree/Orthtree.txt @@ -10,9 +10,9 @@ namespace CGAL { \section Section_Orthtree_Introduction Introduction Quadtrees are tree data structures in which each node encloses a -square section of space, and each internal node has exactly 4 +rectangular section of space, and each internal node has exactly 4 children. Octrees are a similar data structure in 3D in which each -node encloses a cubic section of space, and each internal node has +node encloses a rectangular cuboid section of space, and each internal node has exactly 8 children. We call the generalization of such data structure "orthtrees", as @@ -22,20 +22,32 @@ structures in dimensions 4 and higher. This package provides a general data structure `Orthtree` along with aliases for `Quadtree` and `Octree`. These trees can be constructed -with custom point ranges and split predicates, and iterated on with -various traversal methods. - -\cgalFigureBegin{Orthtree_fig, orthtree.png} -Building an %orthtree in 3D (%octree) from a point cloud. -\cgalFigureEnd - +with custom contents and split predicates, and iterated on with +various traversal methods. Orthants can be orthotopes and not only hypercubes. + +\cgalFigureAnchor{Orthtree_fig } +
+ + + + + + + +
+
+\cgalFigureCaptionBegin{Orthtree_fig} +Top: an %orthtree in 3D (%octree) from a point cloud (top); +Bottom: an %orthtree in 3D (%octree) from the triangle faces of a mesh. +\cgalFigureCaptionEnd \section Section_Orthtree_Building Building -An orthtree is created using a set of points. The points are not -copied: the provided point range is used directly and is rearranged by +A common purpose for an orthtree is to subdivide a collection of points, +and the Orthtree package provides a traits class for this purpose. +The points are not copied: the provided point range is used directly and rearranged by the orthtree. Altering the point range after creating the orthtree -might leave it in an invalid state. The constructor returns a tree +may leave the tree in an invalid state. The constructor returns a tree with a single (root) node that contains all the points. The method [refine()](@ref CGAL::Orthtree::refine) must be called to @@ -43,9 +55,11 @@ subdivide space further. This method uses a split predicate which takes a node as input and returns `true` if this node should be split, `false` otherwise: this enables users to choose on what criterion should the orthtree be refined. Predefined predicates are -provided such as [Maximum_depth](@ref CGAL::Orthtrees::Maximum_depth) or [Maximum_number_of_inliers](@ref CGAL::Orthtrees::Maximum_number_of_inliers). +provided for the point-set orthtree, +including [Maximum_depth](@ref CGAL::Orthtrees::Maximum_depth) +and [Maximum_contained_elements](@ref CGAL::Orthtrees::Maximum_contained_elements). -The simplest way to create an orthtree is using a vector of points. +The simplest way to create a point-set orthtree is from a vector of points. The constructor generally expects a separate point range and map, but the point map defaults to `Identity_property_map` if none is provided. @@ -55,29 +69,27 @@ the existing ones do not match users' needs. \subsection Section_Orthtree_Quadtree Building a Quadtree -The `Orthtree` class may be templated with `Orthtree_traits_2` and thus -behave as a %quadtree. For convenience, the alias `Quadtree` is provided. +The `Orthtree` class may be templated with `Orthtree_traits_point<>` +while specifying a 2d ambient space and thus behave as a %quadtree. +For convenience, the alias `CGAL::Quadtree` is provided. -The following example shows how to create a %quadtree object from a -vector of `Point_2` objects and refine it, which means constructing -the tree's space subdivision itself, using a maximum depth of 10 and a -maximum number of inliers per node (bucket size) of 5. The refinement -is stopped as soon as one of the conditions is violated: if a node has -more inliers than `bucket_size` but is already at `max_depth`, it is -not split. Similarly, a node that is at a depth smaller than -`max_depth` but already has fewer inliers than `bucket_size` is not -split. +The following example shows the construction of %quadtree from a vector of `Point_2` objects. +`quadtree.refine(10, 5)` uses the default split predicate, which +enforces a max-depth and a maximum of elements contained per node ("bucket size"). +Nodes are split if their depth is less than 10, and they contain more than 5 points. \cgalExample{Orthtree/quadtree_build_from_point_vector.cpp} \subsection Section_Orthtree_Point_Vector Building an Octree -The `Orthtree` class may be templated with `Orthtree_traits_3` and thus -behave as an %octree. For convenience, the alias `Octree` is provided. +`Orthtree_traits_point<>` can also be templated with dimension 3 and thus +behave as an %octree. For convenience, the alias `CGAL::Octree` is provided. + +The following example shows how to create an %octree from a vector of `Point_3` objects. +As with the %quadtree example, we use the default split predicate. +In this case the maximum depth is 10, and the bucket size is 1. -The following example shows how to create an %octree from a vector of -`Point_3` objects: \cgalExample{Orthtree/octree_build_from_point_vector.cpp} @@ -87,12 +99,12 @@ Some data structures such as `Point_set_3` require a non-default point map type and object. This example illustrates how to create an octree from a `Point_set_3` loaded from a file. It also shows a more explicit way of setting the split predicate when refining the tree. -An octree is constructed from the point set and its map. -The tree is refined with a maximum depth (deepest node allowed) of 10, -and a bucket size (maximum number of points contained by a single node) of 20. -The tree is then written to the standard output. +An %octree is constructed from the point set and its corresponding map, +and then written to the standard output. The split predicate is manually constructed and passed to the refine method. +In this case, we use a maximum number of contained elements with no corresponding maximum depth. +This means that nodes will continue to be split until none contain more than 10 points. \cgalExample{Orthtree/octree_build_from_point_set.cpp} @@ -108,38 +120,47 @@ at depth 7 can hold 14. \subsection Section_Orthtree_Orthtree_Point_Vector Building an Orthtree -The following example shows how to build an generalized orthtree in dimension 4. -A `std::vector` is manually filled with points. -The vector is used as the point set, -an `Identity_property_map` is automatically set as the orthtree's map type, so a map does not need to be provided. +An orthtree can also be used with an arbitrary number of dimensions. +The `Orthtree_traits_point` template can infer the arbitrary dimension count from the d-dimensional kernel. + +The following example shows how to build a generalized orthtree in dimension 4. +As `std::vector` is manually filled with 4-dimensional points. +The vector is used as the point set, and an `Identity_property_map` is automatically +set as the orthtree's map type, so a map does not need to be provided. \cgalExample{Orthtree/orthtree_build.cpp} -\section Section_Orthtree_Traversal Traversal +\section Section_Orthtree_Properties Properties +The Orthtree uses a mechanism to attach properties to nodes at run-time which follows \ref sectionSurfaceMesh_properties "Surface mesh properties". Each property is identified by a string and its value type and stored in a consecutive block of memory. -\note For simplicity, the rest of the user manual will only use -octrees, but all the presented features also apply to quadtrees and -higher dimension orthtrees. +\section Section_Orthtree_Traversal Traversal %Traversal is the act of navigating among the nodes of the tree. -The `Orthtree` and [Node](@ref CGAL::Orthtree::Node) classes provide a +The `Orthtree` class provides a number of different solutions for traversing the tree. \subsection Section_Orthtree_Manual_Traveral Manual Traversal Because our orthtree is a form of connected acyclic undirected graph, it is possible to navigate between any two nodes. -What that means in practice, is that given a node on the tree, it is possible to +What that means in practice is that given a node on the tree, it is possible to access any other node using the right set of operations. -The `Node` class provides functions that enable the user to access each of its children, as well as its parent (if it exists). +The `Node_index` type provides a handle on a node, and the `orthtree` class provides methods +that enable the user to retrieve the indices of each of its children as well as its parent (if they exist). -The following example demonstrates ways of accessing different nodes of a tree, given a reference to one. +Given any node index `nid`, the `n`th child of that node can be found with `orthtree.child(nid, n)`. +For an octree, values of `n` from 0-7 provide access to the different children. +For non-root nodes, it is possible to access parent nodes using the `orthtree.parent()` accessor. -From the root node, children can be accessed using the subscript operator `CGAL::Orthtree::Node::operator[]()`. -For an octree, values from 0-7 provide access to the different children. +To access grandchildren, it isn't necessary to use nested `orthtree.child()` calls. +Instead, the `orthtree.descendant(node, a, b, ...)` convenience method is provided. +This retrieves the `b`th child of the `a`th child, to any depth. -For non-root nodes, it is possible to access parent nodes using the [parent()](@ref CGAL::Orthtree::Node::parent) accessor. +In most cases, we want to find the descendants of the root node. +For this case, there is another convenience method `orthtree.node(a, b, c, ...)`. +This retrieves the node specified by the descent `a`, `b`, `c`. +It is equivalent to `orthtree.descendant(orthtree.root(), a, b, c, ...)`. -These accessors and operators can be chained to access any node in the tree in a single line of code, as shown in the following example: +The following example demonstrates the use of several of these accessors. \cgalExample{Orthtree/octree_traversal_manual.cpp} @@ -147,9 +168,10 @@ These accessors and operators can be chained to access any node in the tree in a It is often useful to be able to iterate over the nodes of the tree in a particular order. For example, the stream operator `<<` uses a traversal to print out each node. -A few traversals are provided, among them [Preorder_traversal](@ref CGAL::Orthtrees::Preorder_traversal) and [Postorder_traversal](@ref CGAL::Orthtrees::Postorder_traversal). -To traverse a tree in preorder is to visit each parent immediately followed by its children, -whereas in postorder, traversal the children are visited first. +A few traversals are provided, among them [Preorder_traversal](@ref CGAL::Orthtrees::Preorder_traversal) +and [Postorder_traversal](@ref CGAL::Orthtrees::Postorder_traversal). +Traversing a tree in preorder means to visit each parent immediately followed by its children, +whereas in postorder traversal the children are visited first. The following example illustrates how to use the provided traversals. @@ -161,13 +183,15 @@ In this case, we print out the nodes of the tree without indentation instead. \subsection Section_Orthtree_Custom_Traversal Custom Traversal -Users can define their own traversal methods by creating models of the -`OrthtreeTraversal` concept. The following example shows how to define a -custom traversal that only traverses the first branch of the octree: +Users can define their own traversal methods by creating models of the `OrthtreeTraversal` concept. +The custom traversal must provide a method which returns the starting point of the traversal (`first_index()`) +and another method which returns the next node in the sequence (`next_index()`). +`next_index()` returns an empty optional when the end of the traversal is reached. +The following example shows how to define a custom traversal that only traverses the first branch an octree: \cgalExample{Orthtree/octree_traversal_custom.cpp} -\subsection Comparison of Traversals +\subsection Section_Orthtree_Cmp_Trav Comparison of Traversals Figure \cgalFigureRef{Orthtree_traversal_fig} shows in which order nodes are visited depending on the traversal method used. @@ -184,12 +208,11 @@ Once an orthtree is built, its structure can be used to accelerate different tas \subsection Section_Orthtree_Nearest_Neighbor Finding the Nearest Neighbor of a Point -The naive way of finding the nearest neighbor of a point requires finding the distance to every other point. +The naive way of finding the nearest neighbor of a point requires finding the distance to all elements. An orthtree can be used to perform the same task in significantly less time. -For large numbers of points, this can be a large enough difference to outweigh the time spent building the tree. +For large numbers of elements, this can be a large enough difference to outweigh the time spent building the tree. -Note that a kd-tree is expected to outperform the orthtree for this task, -it should be preferred unless features specific to the orthtree are needed. +Note that a kd-tree is expected to outperform the orthtree for this task on points, it should be preferred unless features specific to the orthtree are needed. The following example illustrates how to use an octree to accelerate the search for points close to a location. @@ -200,9 +223,14 @@ Results are put in a vector, and then printed. \cgalExample{Orthtree/octree_find_nearest_neighbor.cpp} +Not all octrees are compatible with nearest neighbor functionality, +as the idea of a nearest neighbor may not make sense for some tree contents. +For the nearest neighbor methods to work, the traits class must implement the +`CollectionPartitioningOrthtreeTraits` concept. + \subsection Section_Orthtree_Grade Grading -An orthtree is graded if the difference of depth between two adjacent +An orthtree is considered "graded" if the difference of depth between two adjacent leaves is at most 1 for every pair of leaves. \cgalFigureBegin{Orthtree_quadree_graded_fig, quadtree_graded.png} @@ -256,15 +284,90 @@ purposes. For nontrivial point counts, the naive approach's calculation time dwarfs that of either the %orthtree or kd-tree. +\section Section_Orthtree_Migration Migrating Code Written Before Release 6.0 + +The orthtree package changed to allow for custom data stored per node in the orthtree. To migrate existing code written before \cgal 6.0 `Orthtree_traits_point` can be used for a point-based orthtrees. The aliases `CGAL::Quadtree` and `CGAL::Octree` have been extended by a boolean template parameter to allow for non-cubic cells, which is the default. The data is passed via the traits class and no longer directly to the orthtree. + +Former code to declare and define an Octree with cubic cells was as follows: +\code{.cpp} +typedef CGAL::Point_set_3 Point_set; +typedef CGAL::Octree Octree; +... +Octree octree(points, points.point_map()); +\endcode + +\cgal 6.0 code with identical behavior is now: +\code{.cpp} +typedef CGAL::Point_set_3 Point_set; +typedef CGAL::Octree Octree; +... +Octree octree(points, points.point_map()); +\endcode + +The node class does not exist anymore and has been replaced by the lightweight type `Node_index`. All information formerly contained in the node class is now accessible via the `Orthtree` interface. + +Former node access was as follows: +\code{.cpp} +Orthtree::Node root = orthtree.root(); +Orthtree::Node child = root[0]; +bool is_leaf = child.is_leaf(); + +for (Orthtree::Node::const_iterator it = child.begin(); it != child.end(); it++) { + const Orthtree::Point &p = *it; + ... +} +\endcode + +\cgal 6.0 node access is now: +\code{.cpp} +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +Point_set points; +... +Orthtree::Node_index root = orthtree.root(); +Orthtree::Node_index child = orthtree.child(root, 0); +bool is_leaf = orthtree.is_leaf(child); + +for (Octree::Traits::Node_data_element& e : octree.data(child)) { + // Using a pointmap is necessary when using a Point_set_3. + Point& p = points.point(e); + + // If the orthtree is used with a std::vector, Node_data_element is identical to Point. + // Point& p = e; + ... +} +\endcode + +The nearest neighbor search behaves as before, however, the output iterator will be required to store `CollectionPartitioningOrthtreeTraits::Node_data_element`. This may be the actual point or, e.g., in case a `Point_set_3` is used, an index which requires the use of a point map to retrieve the point. + +The provided traversals, i.e., `CGAL::Orthtrees::Leaves_traversal`, `CGAL::Orthtrees::Preorder_traversal`, `CGAL::Orthtrees::Postorder_traversal`, require the orthtree as template parameter now. + +Former traversal use was as follows: +\code{.cpp} +for (Orthtree::Node node : orthtree.traverse()) + ... +\endcode + +\cgal 6.0 traversal use is now: +\code{.cpp} +for (Orthtree::Node_index node : orthtree.traverse>()) + ... +\endcode + \section Section_Orthtree_History History A prototype code was implemented by Pierre Alliez and improved by Tong Zhao and Cédric Portaneri. From this prototype code, the package was -developed by Jackson Campolatarro as part of the Google Summer of Code -2020. Simon Giraudot, supervisor of the GSoC internship, completed and +developed by Jackson Campolattaro as part of the Google Summer of Code 2020. +Simon Giraudot, supervisor of the GSoC internship, completed and finalized the package for integration in CGAL 5.3. Pierre Alliez provided kind help and advice all the way through. +Starting with CGAL 6.0 an API improvement of the Orthtree package was released. +Most notably, the orthtree nodes do not need to store anything. Data to be stored +by the node are managed through a mechanism of dynamic properties. +This improvement was done by Jackson Campolattaro thanks to a funding provided by +INRIA, together with GeometryFactory. */ -} +} \ No newline at end of file diff --git a/Orthtree/doc/Orthtree/PackageDescription.txt b/Orthtree/doc/Orthtree/PackageDescription.txt index 22b02350771a..98825d8e237d 100644 --- a/Orthtree/doc/Orthtree/PackageDescription.txt +++ b/Orthtree/doc/Orthtree/PackageDescription.txt @@ -1,20 +1,20 @@ -/// \defgroup PkgOrthtreeRef Quadtree\, Octree and Orthtree Reference +/// \defgroup PkgOrthtreeRef Quadtrees, Octrees and Orthtrees Reference +Quadtree, Octree and Orthtree Reference -/// \defgroup PkgOrthtreeClasses Classes +/// \defgroup PkgOrthtreeConcepts Concepts /// \ingroup PkgOrthtreeRef /// \defgroup PkgOrthtreeTraits Traits -/// \ingroup PkgOrthtreeClasses +/// \ingroup PkgOrthtreeRef /// \defgroup PkgOrthtreeSplitPredicates Split Predicates -/// \ingroup PkgOrthtreeClasses +/// \ingroup PkgOrthtreeRef /// \defgroup PkgOrthtreeTraversal Traversal -/// \ingroup PkgOrthtreeClasses - -/// \defgroup PkgOrthtreeConcepts Concepts /// \ingroup PkgOrthtreeRef + + /*! \addtogroup PkgOrthtreeRef @@ -39,22 +39,25 @@ \cgalCRPSection{Concepts} - `OrthtreeTraits` +- `OrthtreeTraitsWithData` - `OrthtreeTraversal` +- `CollectionPartitioningOrthtreeTraits` \cgalCRPSection{Classes} - `CGAL::Quadtree` - `CGAL::Octree` -- `CGAL::Orthtree` +- `CGAL::Orthtree` \cgalCRPSection{Traits} -- `CGAL::Orthtree_traits_2` -- `CGAL::Orthtree_traits_3` -- `CGAL::Orthtree_traits_d` +- `CGAL::Orthtree_traits` +- `CGAL::Orthtree_traits_point` +- `CGAL::Orthtree_traits_base` +- `CGAL::Orthtree_traits_face_graph` \cgalCRPSection{Split Predicates} -- `CGAL::Orthtrees::Maximum_number_of_inliers` +- `CGAL::Orthtrees::Maximum_contained_elements` - `CGAL::Orthtrees::Maximum_depth` -- `CGAL::Orthtrees::Maximum_depth_and_maximum_number_of_inliers` +- `CGAL::Orthtrees::Maximum_depth_and_maximum_contained_elements` \cgalCRPSection{%Traversal} - `CGAL::Orthtrees::Preorder_traversal` diff --git a/Orthtree/doc/Orthtree/dependencies b/Orthtree/doc/Orthtree/dependencies index 1e72914f7008..703f46d89264 100644 --- a/Orthtree/doc/Orthtree/dependencies +++ b/Orthtree/doc/Orthtree/dependencies @@ -8,3 +8,4 @@ Property_map STL_Extension Spatial_searching Stream_support +Surface_mesh diff --git a/Orthtree/doc/Orthtree/examples.txt b/Orthtree/doc/Orthtree/examples.txt index ff327b8272fd..7d6ccf98a35f 100644 --- a/Orthtree/doc/Orthtree/examples.txt +++ b/Orthtree/doc/Orthtree/examples.txt @@ -3,6 +3,7 @@ \example Orthtree/octree_build_from_point_vector.cpp \example Orthtree/octree_build_with_custom_split.cpp \example Orthtree/octree_find_nearest_neighbor.cpp +\example Orthtree/octree_surface_mesh.cpp \example Orthtree/octree_traversal_manual.cpp \example Orthtree/octree_traversal_preorder.cpp \example Orthtree/octree_traversal_custom.cpp diff --git a/Orthtree/doc/Orthtree/fig/orthtree_mesh.png b/Orthtree/doc/Orthtree/fig/orthtree_mesh.png new file mode 100644 index 000000000000..12864286cd0a Binary files /dev/null and b/Orthtree/doc/Orthtree/fig/orthtree_mesh.png differ diff --git a/Orthtree/examples/Orthtree/CMakeLists.txt b/Orthtree/examples/Orthtree/CMakeLists.txt index 432d99b21c48..ee2c7502a8da 100644 --- a/Orthtree/examples/Orthtree/CMakeLists.txt +++ b/Orthtree/examples/Orthtree/CMakeLists.txt @@ -15,6 +15,8 @@ create_single_source_cgal_program("octree_traversal_manual.cpp") create_single_source_cgal_program("octree_traversal_preorder.cpp") create_single_source_cgal_program("octree_grade.cpp") create_single_source_cgal_program("quadtree_build_from_point_vector.cpp") +create_single_source_cgal_program("octree_surface_mesh.cpp") +create_single_source_cgal_program("quadtree_build_manually.cpp") find_package(Eigen3 3.1.91 QUIET) #(requires 3.1.91 or greater) include(CGAL_Eigen_support) diff --git a/Orthtree/examples/Orthtree/octree_build_from_point_set.cpp b/Orthtree/examples/Orthtree/octree_build_from_point_set.cpp index 4ef69692131e..c4226b0dbe8a 100644 --- a/Orthtree/examples/Orthtree/octree_build_from_point_set.cpp +++ b/Orthtree/examples/Orthtree/octree_build_from_point_set.cpp @@ -7,12 +7,11 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; int main(int argc, char **argv) { @@ -33,10 +32,10 @@ int main(int argc, char **argv) { Octree octree(points, points.point_map()); // Build the octree with a small bucket size, using a more verbose method - octree.refine(CGAL::Orthtrees::Maximum_depth_and_maximum_number_of_inliers(5, 10)); + octree.refine(CGAL::Orthtrees::Maximum_contained_elements(10)); // Print out the tree - std::cout << octree; + std::cout << octree << std::endl; return EXIT_SUCCESS; } diff --git a/Orthtree/examples/Orthtree/octree_build_from_point_vector.cpp b/Orthtree/examples/Orthtree/octree_build_from_point_vector.cpp index 71f1393df32b..1081e155df30 100644 --- a/Orthtree/examples/Orthtree/octree_build_from_point_vector.cpp +++ b/Orthtree/examples/Orthtree/octree_build_from_point_vector.cpp @@ -4,11 +4,10 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef std::list Point_vector; - -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_vector = std::vector; +using Octree = CGAL::Octree; int main() { @@ -27,7 +26,7 @@ int main() { Octree octree(points); // Build the octree - octree.refine(10, 20); + octree.refine(10, 1); // Print out the tree std::cout << octree; diff --git a/Orthtree/examples/Orthtree/octree_build_with_custom_split.cpp b/Orthtree/examples/Orthtree/octree_build_with_custom_split.cpp index 3202c6c9ade1..92c7c777a1d0 100644 --- a/Orthtree/examples/Orthtree/octree_build_with_custom_split.cpp +++ b/Orthtree/examples/Orthtree/octree_build_with_custom_split.cpp @@ -7,25 +7,26 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; // Split Predicate -// The predicate is a functor which returns a boolean value, whether a node needs to be split or not +// The predicate is a functor which returns a Boolean value, whether a node needs to be split or not struct Split_by_ratio { std::size_t ratio; - Split_by_ratio(std::size_t ratio) : ratio(ratio) {} + explicit Split_by_ratio(std::size_t ratio) : ratio(ratio) {} - template - bool operator()(const Node &n) const { - return n.size() > (ratio * n.depth()); + template + bool operator()(Node_index i, const Tree &tree) const { + std::size_t num_points = tree.data(i).size(); + std::size_t depth = tree.depth(i); + return num_points > (ratio * depth); } }; diff --git a/Orthtree/examples/Orthtree/octree_find_nearest_neighbor.cpp b/Orthtree/examples/Orthtree/octree_find_nearest_neighbor.cpp index e94fca1ee206..2d687aa42159 100644 --- a/Orthtree/examples/Orthtree/octree_find_nearest_neighbor.cpp +++ b/Orthtree/examples/Orthtree/octree_find_nearest_neighbor.cpp @@ -9,14 +9,13 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; -typedef CGAL::Octree Octree; - -int main(int argc, char **argv) { +int main(int argc, char** argv) { // Point set will be used to hold our points Point_set points; @@ -39,21 +38,42 @@ int main(int argc, char **argv) { // Find the nearest points to a few locations std::vector points_to_find = { - {0, 0, 0}, - {1, 1, 1}, - {-1, -1, -1}, - {-0.46026, -0.25353, 0.32051}, - {-0.460261, -0.253533, 0.320513} + {0, 0, 0}, + {1, 1, 1}, + {-1, -1, -1}, + {-0.46026, -0.25353, 0.32051}, + {-0.460261, -0.253533, 0.320513} }; + for (const Point& p : points_to_find) - octree.nearest_neighbors - (p, 1, // k=1 to find the single closest point - boost::make_function_output_iterator - ([&](const Point& nearest) + octree.nearest_k_neighbors + (p, 1, // k=1 to find the single closest point + boost::make_function_output_iterator + ([&](const Point_set::Index& nearest) { std::cout << "the nearest point to (" << p << - ") is (" << nearest << ")" << std::endl; + ") is (" << points.point(nearest) << ")" << std::endl; })); + typename Octree::Sphere s(points_to_find[0], 0.0375); + std::cout << std::endl << "Closest points within the sphere around " << s.center() << " with squared radius of " << s.squared_radius() << ":" << std::endl; + octree.neighbors_within_radius + (s, + boost::make_function_output_iterator + ([&](const Point_set::Index& nearest) + { + std::cout << points.point(nearest) << " dist: " << (Point(0, 0, 0) - points.point(nearest)).squared_length() << std::endl; + })); + + std::size_t k = 3; + std::cout << std::endl << "The up to " << k << " closest points to(" << s.center() << ") within a squared radius of " << s.squared_radius() << " are:" << std::endl; + octree.nearest_k_neighbors_within_radius + (s, k, + boost::make_function_output_iterator + ([&](const Point_set::Index& nearest) + { + std::cout << "(" << points.point(nearest) << " dist: " << (Point(0, 0, 0) - points.point(nearest)).squared_length() << std::endl; + })); + return EXIT_SUCCESS; } diff --git a/Orthtree/examples/Orthtree/octree_grade.cpp b/Orthtree/examples/Orthtree/octree_grade.cpp index b00011045c68..5f213ff1233e 100644 --- a/Orthtree/examples/Orthtree/octree_grade.cpp +++ b/Orthtree/examples/Orthtree/octree_grade.cpp @@ -4,10 +4,10 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef std::vector Point_vector; -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_vector = std::vector; +using Octree = CGAL::Octree; int main() { diff --git a/Orthtree/examples/Orthtree/octree_surface_mesh.cpp b/Orthtree/examples/Orthtree/octree_surface_mesh.cpp new file mode 100644 index 000000000000..676eacae2739 --- /dev/null +++ b/Orthtree/examples/Orthtree/octree_surface_mesh.cpp @@ -0,0 +1,68 @@ +#include +#include + +#include +#include + +#include + +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Mesh = CGAL::Surface_mesh; + +using OTraits = CGAL::Orthtree_traits_face_graph>; +using Octree = CGAL::Orthtree; + +void dump_as_polylines(const Octree& ot) +{ + std::ofstream out("octree.polylines.txt"); + for (Octree::Node_index node : ot.traverse(CGAL::Orthtrees::Leaves_traversal(ot))) + { + if (!ot.is_leaf(node)) + continue; + auto bb = ot.bbox(node); + out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin() + << " " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin() << "\n"; + out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin() + << " " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin() << "\n"; + out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin() + << " " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n"; + out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin() + << " " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() << "\n"; +// + out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmin() + << " " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmin() + << " " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmin() + << " " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmax() << " " << bb.ymax() << " " << bb.zmin() + << " " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n"; +// + out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax() + << " " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmin() << " " << bb.ymin() << " " << bb.zmax() + << " " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmax() << " " << bb.ymin() << " " << bb.zmax() + << " " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n"; + out << "2 " << bb.xmin() << " " << bb.ymax() << " " << bb.zmax() + << " " << bb.xmax() << " " << bb.ymax() << " " << bb.zmax() << "\n"; + } +} + +int main(int argc, char** argv) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off"); + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh)) + { + std::cerr << "Error: cannot read file" << std::endl; + return EXIT_FAILURE; + } + + Octree tree(mesh, mesh.points()); + OTraits::Split_predicate_node_min_extent sp(0.01); + tree.refine(sp); + + dump_as_polylines(tree); +} diff --git a/Orthtree/examples/Orthtree/octree_traversal_custom.cpp b/Orthtree/examples/Orthtree/octree_traversal_custom.cpp index 3347b54ff996..a9a3ccba0728 100644 --- a/Orthtree/examples/Orthtree/octree_traversal_custom.cpp +++ b/Orthtree/examples/Orthtree/octree_traversal_custom.cpp @@ -3,34 +3,42 @@ #include #include +#include #include #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; -typedef Octree::Node Node; - -struct First_branch_traversal -{ - Node first (Node root) const - { - return root; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; + +template +struct First_branch_traversal { + + using Node_index = typename Tree::Node_index; + + const Tree& m_orthtree; + + explicit First_branch_traversal(const Tree& orthtree) : m_orthtree(orthtree) {} + + Node_index first_index() const { + return m_orthtree.root(); } - Node next (Node n) const - { - if (n.is_leaf()) - return Node(); - return n[0]; + std::optional next_index(Node_index n) const { + + // Stop when we reach the base of the tree + if (m_orthtree.is_leaf(n)) + return {}; + + // Always descend the first child + return m_orthtree.child(n, 0); } }; -int main(int argc, char **argv) { +int main(int argc, char** argv) { // Point set will be used to hold our points Point_set points; @@ -52,8 +60,9 @@ int main(int argc, char **argv) { octree.refine(); // Print out the first branch using custom traversal - for (auto &node : octree.traverse()) - std::cout << node << std::endl; + for (auto node: octree.traverse>()) { + std::cout << octree.to_string(node) << std::endl; + } return EXIT_SUCCESS; } diff --git a/Orthtree/examples/Orthtree/octree_traversal_manual.cpp b/Orthtree/examples/Orthtree/octree_traversal_manual.cpp index dcbe8d0f15c0..45c51fa3797e 100644 --- a/Orthtree/examples/Orthtree/octree_traversal_manual.cpp +++ b/Orthtree/examples/Orthtree/octree_traversal_manual.cpp @@ -7,14 +7,13 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; -typedef CGAL::Octree Octree; - -int main(int argc, char **argv) { +int main(int argc, char** argv) { // Point set will be used to hold our points Point_set points; @@ -39,29 +38,30 @@ int main(int argc, char **argv) { std::cout << "Navigation relative to the root node" << std::endl; std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; std::cout << "the root node: " << std::endl; - std::cout << octree.root() << std::endl; + std::cout << octree.to_string(octree.root()) << std::endl; std::cout << "the first child of the root node: " << std::endl; - std::cout << octree.root()[0] << std::endl; + std::cout << octree.to_string(octree.child(octree.root(), 0)) << std::endl; std::cout << "the fifth child: " << std::endl; - std::cout << octree.root()[4] << std::endl; - std::cout << "the fifth child, accessed without the root keyword: " << std::endl; - std::cout << octree[4] << std::endl; + std::cout << octree.to_string(octree.child(octree.root(), 4)) << std::endl; + std::cout << "the fifth child, accessed without going through root: " << std::endl; + std::cout << octree.to_string(octree.node(4)) << std::endl; std::cout << "the second child of the fourth child: " << std::endl; - std::cout << octree.root()[4][1] << std::endl; - std::cout << "the second child of the fourth child, accessed without the root keyword: " << std::endl; - std::cout << octree[4][1] << std::endl; + std::cout << octree.to_string(octree.child(octree.child(octree.root(), 4), 1)) << std::endl; + std::cout << "the second child of the fourth child, accessed without going through root: " << std::endl; + std::cout << octree.to_string(octree.node(4, 1)) << std::endl; std::cout << std::endl; // Retrieve one of the deeper children - Octree::Node cur = octree[3][2]; + Octree::Node_index cur = octree.node(4, 3); std::cout << "Navigation relative to a child node" << std::endl; std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; std::cout << "the third child of the fourth child: " << std::endl; - std::cout << cur << std::endl; + std::cout << octree.to_string(cur) << std::endl; std::cout << "the third child: " << std::endl; - std::cout << cur.parent() << std::endl; + std::cout << octree.to_string(octree.parent(cur)) << std::endl; std::cout << "the next sibling of the third child of the fourth child: " << std::endl; - std::cout << cur.parent()[cur.local_coordinates().to_ulong() + 1] << std::endl; + std::cout << octree.to_string(octree.child(octree.parent(cur), octree.local_coordinates(cur).to_ulong() + 1)) + << std::endl; return EXIT_SUCCESS; } diff --git a/Orthtree/examples/Orthtree/octree_traversal_preorder.cpp b/Orthtree/examples/Orthtree/octree_traversal_preorder.cpp index 72587767a607..a359f6217eea 100644 --- a/Orthtree/examples/Orthtree/octree_traversal_preorder.cpp +++ b/Orthtree/examples/Orthtree/octree_traversal_preorder.cpp @@ -7,13 +7,12 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef Point_set::Point_map Point_map; - -typedef CGAL::Octree Octree; -typedef CGAL::Orthtrees::Preorder_traversal Preorder_traversal; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Point_map = Point_set::Point_map; +using Octree = CGAL::Octree; +using Preorder_traversal = CGAL::Orthtrees::Preorder_traversal; int main(int argc, char **argv) { @@ -37,9 +36,8 @@ int main(int argc, char **argv) { octree.refine(); // Print out the octree using preorder traversal - for (Octree::Node node : octree.traverse()) { - - std::cout << node << std::endl; + for (auto node : octree.traverse()) { + std::cout << octree.to_string(node) << std::endl; } return EXIT_SUCCESS; diff --git a/Orthtree/examples/Orthtree/orthtree_build.cpp b/Orthtree/examples/Orthtree/orthtree_build.cpp index f553eccc695c..56419faf0271 100644 --- a/Orthtree/examples/Orthtree/orthtree_build.cpp +++ b/Orthtree/examples/Orthtree/orthtree_build.cpp @@ -1,27 +1,27 @@ #include #include +#include #include -#include +#include #include // Type Declarations -typedef CGAL::Dimension_tag<4> Dimension; -typedef CGAL::Epick_d Kernel; -typedef Kernel::Point_d Point_d; -typedef std::vector Point_vector; - -typedef CGAL::Orthtree_traits_d Traits; -typedef CGAL::Orthtree Orthtree; +const int dimension = 4; +using Kernel = CGAL::Epick_d >; +using Point_d = Kernel::Point_d; +using Point_vector = std::vector; +using Traits = CGAL::Orthtree_traits_point; +using Orthtree = CGAL::Orthtree; int main() { CGAL::Random r; Point_vector points_dd; - for (std::size_t i = 0; i < 5; ++ i) + for (std::size_t i = 0; i < 20; ++ i) { - std::array init; + std::array init{}; for (double& v : init) v = r.get_double(-1., 1.); points_dd.emplace_back (init.begin(), init.end()); @@ -30,5 +30,7 @@ int main() Orthtree orthtree(points_dd); orthtree.refine(10, 5); + std::cout << orthtree << std::endl; + return EXIT_SUCCESS; } diff --git a/Orthtree/examples/Orthtree/quadtree_build_from_point_vector.cpp b/Orthtree/examples/Orthtree/quadtree_build_from_point_vector.cpp index 345652688572..bd40ace04232 100644 --- a/Orthtree/examples/Orthtree/quadtree_build_from_point_vector.cpp +++ b/Orthtree/examples/Orthtree/quadtree_build_from_point_vector.cpp @@ -5,11 +5,10 @@ #include // Type Declarations -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_2 Point_2; -typedef std::vector Point_vector; - -typedef CGAL::Quadtree Quadtree; +using Kernel = CGAL::Simple_cartesian; +using Point_2 = Kernel::Point_2; +using Point_vector = std::vector; +using Quadtree = CGAL::Quadtree; int main() { diff --git a/Orthtree/examples/Orthtree/quadtree_build_manually.cpp b/Orthtree/examples/Orthtree/quadtree_build_manually.cpp new file mode 100644 index 000000000000..2e38f406c384 --- /dev/null +++ b/Orthtree/examples/Orthtree/quadtree_build_manually.cpp @@ -0,0 +1,59 @@ +#include + +#include +#include + +#include +#include + +using Kernel = CGAL::Simple_cartesian; + +namespace CGAL { + +struct empty_type { +}; + +template +struct Orthtree_traits_empty : public Orthtree_traits_base { + + using Self = Orthtree_traits_empty; + using Tree = Orthtree; + + using Node_data = std::array; + + Orthtree_traits_empty(typename Self::Bbox_d bbox) : m_bbox(bbox) {}; + + auto construct_root_node_bbox_object() const { + return [&]() -> typename Self::Bbox_d { return m_bbox; }; + } + + auto construct_root_node_contents_object() const { return [&]() -> Node_data { return {}; }; } + + auto distribute_node_contents_object() { + return [&](typename Tree::Node_index /* n */, Tree& /* tree */, const typename Self::Point_d& /* center */) -> void {}; + } + +private: + + typename Self::Bbox_d m_bbox; + +}; +} + +using EmptyQuadtree = CGAL::Orthtree>; + +int main() { + + // Build an empty quadtree which covers the domain (-1, -1) to (1, 1) + EmptyQuadtree quadtree{{{-1, -1, 1, 1}}}; + + // Split several nodes of the tree + quadtree.split(quadtree.root()); + quadtree.split(quadtree.node(0)); + quadtree.split(quadtree.node(3)); + quadtree.split(quadtree.node(3, 0)); + + std::cout << quadtree.bbox(quadtree.root()) << std::endl; + std::cout << quadtree << std::endl; + return EXIT_SUCCESS; +} diff --git a/Orthtree/include/CGAL/Octree.h b/Orthtree/include/CGAL/Octree.h index bd7e50b3b881..e59a19be46d5 100644 --- a/Orthtree/include/CGAL/Octree.h +++ b/Orthtree/include/CGAL/Octree.h @@ -15,34 +15,27 @@ #include #include -#include +#include namespace CGAL { /*! - \ingroup PkgOrthtreeClasses + \ingroup PkgOrthtreeRef - \brief Alias that specializes the `Orthtree` class to a 3D octree. + \brief Alias that specializes the `Orthtree` class to a 3D octree storing 3D points. - These two types are exactly equivalent: - - `Octree` - - `Orthtree, PointRange, PointMap>`. - - \warning This is a not a real class but an alias, please refer to - the documentation of `Orthtree`. - - \tparam GeomTraits must be a model of `Kernel` - \tparam PointRange_ must be a model of range whose value type is the key type of `PointMap` - \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_3` + \tparam GeomTraits a model of `Kernel` + \tparam PointRange a model of `Range` whose value type is the key type of `PointMap` + \tparam PointMap a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_3` + \tparam cubic_nodes Boolean to enforce cubic nodes */ -template ::value_type> > -#ifdef DOXYGEN_RUNNING -class Octree; -#else -using Octree = Orthtree, PointRange, PointMap>; -#endif +template < + typename GeomTraits, + typename PointRange, + typename PointMap = Identity_property_map::value_type>, + bool cubic_nodes = false +> +using Octree = Orthtree>; } // namespace CGAL diff --git a/Orthtree/include/CGAL/Orthtree.h b/Orthtree/include/CGAL/Orthtree.h index afc958b97547..0cba3e781ba7 100644 --- a/Orthtree/include/CGAL/Orthtree.h +++ b/Orthtree/include/CGAL/Orthtree.h @@ -20,10 +20,12 @@ #include #include +#include +#include +#include #include #include #include -#include #include #include @@ -39,123 +41,193 @@ #include #include #include +#include + +#include namespace CGAL { +namespace Orthtree_impl { + +BOOST_MPL_HAS_XXX_TRAIT_DEF(Node_data) +BOOST_MPL_HAS_XXX_TRAIT_DEF(Squared_distance_of_element) + +template +struct Node_data_wrapper; + +template +struct Node_data_wrapper +{ + using Node_index = typename GT::Node_index; + using Node_data = typename GT::Node_data; + typename CGAL::Properties::Experimental::Property_container::template Array& m_node_contents; + + template + Node_data_wrapper(Property_container& node_properties) + : m_node_contents(node_properties.template get_or_add_property("contents").first) + {} + + const Node_data& operator[](Node_index n) const + { + return m_node_contents[n]; + } + + Node_data& operator[](Node_index n) + { + return m_node_contents[n]; + } +}; + +template +struct Node_data_wrapper +{ + using Node_index = typename GT::Node_index; + using Node_data = void*; + + template + Node_data_wrapper(Property_container&) {} + + void* operator[](Node_index) const + { + return nullptr; + } +}; + +} // end of Orthtree_impl namespace + /*! - \ingroup PkgOrthtreeClasses + \ingroup PkgOrthtreeRef - \brief A data structure using an axis-aligned hybercubic - decomposition of dD space for efficient point access and - computations. + \brief A data structure using an axis-aligned hyperrectangle + decomposition of dD space for efficient access and + computation. - \details It builds a hierarchy of nodes which subdivide the space - based on a collection of points. Each node represents an - axis-aligned hypercubic region of space. A node contains the range - of points that are present in the region it defines, and it may + \details It builds a hierarchy of nodes which subdivides the space. + Each node represents an axis-aligned hyperrectangle region of space. + The contents of nodes depend on the traits class, non-leaf nodes also contain \f$2^{dim}\f$ other nodes which further subdivide the region. \sa `CGAL::Quadtree` \sa `CGAL::Octree` - \tparam Traits_ must be a model of `OrthtreeTraits` - \tparam PointRange_ must be a model of range whose value type is the key type of `PointMap_` - \tparam PointMap_ must be a model of `ReadablePropertyMap` whose value type is `Traits_::Point_d` + \tparam GeomTraits must be a model of `OrthtreeTraits` or `OrthtreeTraitsWithData`. */ -template > -class Orthtree -{ - +template +class Orthtree { public: - /// \name Template Types /// @{ - typedef Traits_ Traits; ///< Geometry traits - typedef PointRange_ PointRange; ///< Point range - typedef PointMap_ PointMap; ///< Point map + using Traits = GeomTraits; ///< Geometry traits /// @} /// \name Traits Types /// @{ - typedef typename Traits::Dimension Dimension; ///< Dimension of the tree - typedef typename Traits::FT FT; ///< Number type. - typedef typename Traits::Point_d Point; ///< Point type. - typedef typename Traits::Bbox_d Bbox; ///< Bounding box type. - typedef typename Traits::Sphere_d Sphere; ///< Sphere type. +#ifndef DOXYGEN_RUNNING + static inline constexpr bool has_data = Orthtree_impl::has_Node_data::value; + static inline constexpr bool supports_neighbor_search = true;// Orthtree_impl::has_Squared_distance_of_element::value; +#else + static inline constexpr bool has_data = bool_value; ///< `true` if `GeomTraits` is a model of `OrthtreeTraitsWithData` and `false` otherwise. + static inline constexpr bool supports_neighbor_search = bool_value; ///< `true` if `GeomTraits` is a model of `CollectionPartitioningOrthtreeTraits` and `false` otherwise. +#endif + static constexpr int dimension = Traits::dimension; ///< Dimension of the tree + using Kernel = typename Traits::Kernel; ///< Kernel type. + using FT = typename Traits::FT; ///< Number type. + using Point = typename Traits::Point_d; ///< Point type. + using Bbox = typename Traits::Bbox_d; ///< Bounding box type. + using Sphere = typename Traits::Sphere_d; ///< Sphere type. + using Adjacency = typename Traits::Adjacency; ///< Adjacency type. + + using Node_index = typename Traits::Node_index; ///< Index of a given node in the tree; the root always has index 0. +#ifndef DOXYGEN_RUNNING + using Node_data = typename Orthtree_impl::Node_data_wrapper::Node_data; +#else + using Node_data = std::conditional_t; +#endif - /// \cond SKIP_IN_MANUAL - typedef typename Traits::Array Array; ///< Array type. - typedef typename Traits::Construct_point_d_from_array - Construct_point_d_from_array; - typedef typename Traits::Construct_bbox_d - Construct_bbox_d; - /// \endcond /// @} /// \name Public Types /// @{ /*! - * \brief Self typedef for convenience. + * \brief Self alias for convenience. */ - typedef Orthtree Self; + using Self = Orthtree; /*! * \brief Degree of the tree (number of children of non-leaf nodes). */ - typedef Dimension_tag<(2 << (Dimension::value-1))> Degree; + static constexpr int degree = (2 << (dimension - 1)); + + /*! + \brief Set of bits representing this node's relationship to its parent. + + Equivalent to an array of Booleans, where index[0] is whether `x` + is greater, index[1] is whether `y` is greater, index[2] is whether + `z` is greater, and so on for higher dimensions if needed. + Used to represent a node's relationship to the center of its parent. + */ + using Local_coordinates = std::bitset; /*! - * \brief The Sub-tree / Orthant type. + \brief Coordinates representing this node's relationship + with the rest of the tree. + + Each value `(x, y, z, ...)` of global coordinates is calculated by doubling + the parent's global coordinates and adding the local coordinates. */ - class Node; + using Global_coordinates = std::array; /*! * \brief A predicate that determines whether a node must be split when refining a tree. */ - typedef std::function Split_predicate; + using Split_predicate = std::function; /*! - * \brief A model of `ConstRange` whose value type is `Node`. + * \brief A model of `ForwardRange` whose value type is `Node_index`. */ #ifdef DOXYGEN_RUNNING - typedef unspecified_type Node_range; + using Node_index_range = unspecified_type; #else - typedef boost::iterator_range > Node_range; + using Node_index_range = boost::iterator_range>; #endif - /// \cond SKIP_IN_MANUAL - /*! - * \brief A function that determines the next node in a traversal given the current one. + * \brief A model of `LvaluePropertyMap` with `Node_index` as key type and `T` as value type. */ - typedef std::function Node_traversal_method_const; - - /// \endcond - - /// \cond SKIP_IN_MANUAL - typedef typename PointRange::iterator Range_iterator; - typedef typename std::iterator_traits::value_type Range_type; - typedef Orthtrees::internal::Cartesian_ranges Cartesian_ranges; - /// \endcond +#ifdef DOXYGEN_RUNNING + template + using Property_map = unspecified_type; +#else + template + using Property_map = Properties::Experimental::Property_array_handle; +#endif /// @} private: // data members : - Traits m_traits; /* the tree traits */ - PointRange& m_range; /* input point range */ - PointMap m_point_map; /* property map: `value_type of InputIterator` -> `Point` (Position) */ + using Cartesian_ranges = Orthtrees::internal::Cartesian_ranges; + using Node_property_container = Properties::Experimental::Property_container; - Node m_root; /* root node of the orthtree */ + template + using Property_array = typename Properties::Experimental::Property_container::template Array; - Point m_bbox_min; /* input bounding box min value */ + Traits m_traits; /* the tree traits */ + + Node_property_container m_node_properties; + Orthtree_impl::Node_data_wrapper m_node_contents; + Property_array& m_node_depths; + Property_array& m_node_coordinates; + Property_array>& m_node_parents; + Property_array>& m_node_children; - std::vector m_side_per_depth; /* side length per node's depth */ + using Bbox_dimensions = std::array; + Bbox m_bbox; + std::vector m_side_per_depth; /* precomputed (potentially approximated) side length per node's depth */ - Cartesian_ranges cartesian_range; /* a helper to easily iterator on coordinates of points */ + Cartesian_ranges cartesian_range; /* a helper to easily iterate over coordinates of points */ public: @@ -163,141 +235,88 @@ class Orthtree /// @{ /*! - \brief creates an orthtree from a collection of points. - - The constructed orthtree has a root node, with no children, that - contains the points passed. That root node has a bounding box that - encloses all of the points passed, with padding according to the - `enlarge_ratio`. + \brief constructs an orthtree for a traits instance. - That root node is built as a `n`-cube (square in 2D, cube in 3D, - etc.) whose edge size is the longest bounding box edge multiplied - by `enlarge_ratio`. Using an `enlarge_ratio>1.0` prevents some - points from being exactly on the border of some cells, which can - lead to over-refinement. + The constructed orthtree has a root node with no children, + containing the contents determined by `Construct_root_node_contents` from the traits class. + That root node has a bounding box determined by `Construct_root_node_bbox` from the traits class, + which typically encloses its contents. This single-node orthtree is valid and compatible - with all Orthtree functionality, but any performance benefits are + with all orthtree functionality, but any performance benefits are unlikely to be realized until `refine()` is called. - \warning The input point range is not copied. It is used directly - and is rearranged by the `Orthtree`. Altering the point range - after creating the orthtree might leave it in an invalid state. - - \param point_range input point range. - \param point_map property map to access the input points. - \param enlarge_ratio ratio to which the bounding box should be enlarged. \param traits the traits object. */ - Orthtree(PointRange& point_range, - PointMap point_map = PointMap(), - const FT enlarge_ratio = 1.2, - Traits traits = Traits()) - : m_traits (traits) - , m_range (point_range) - , m_point_map (point_map) - , m_root(Node(), 0) - { - Array bbox_min; - Array bbox_max; + explicit Orthtree(Traits traits) : + m_traits(traits), + m_node_contents(m_node_properties), + m_node_depths(m_node_properties.template get_or_add_property("depths", 0).first), + m_node_coordinates(m_node_properties.template get_or_add_property("coordinates").first), + m_node_parents(m_node_properties.template get_or_add_property>("parents").first), + m_node_children(m_node_properties.template get_or_add_property>("children").first) { - // init bbox with first values found - { - const Point& p = get (m_point_map, *(point_range.begin())); - std::size_t i = 0; - for (const FT& x : cartesian_range(p)) - { - bbox_min[i] = x; - bbox_max[i] = x; - ++ i; - } - } + m_node_properties.emplace(); - for (const Range_type& r : point_range) - { - const Point& p = get (m_point_map, r); - std::size_t i = 0; - for (const FT& x : cartesian_range(p)) - { - bbox_min[i] = (std::min)(x, bbox_min[i]); - bbox_max[i] = (std::max)(x, bbox_max[i]); - ++ i; - } - } + // init bbox with first values found + m_bbox = m_traits.construct_root_node_bbox_object()(); - Array bbox_centroid; - FT max_length = FT(0); - for (std::size_t i = 0 ; i < Dimension::value; ++ i) - { - bbox_centroid[i] = (bbox_min[i] + bbox_max[i]) / FT(2); - max_length = (std::max)(max_length, bbox_max[i] - bbox_min[i]); - } - max_length *= enlarge_ratio / FT(2); + // Determine dimensions of the root bbox - for (std::size_t i = 0 ; i < Dimension::value; ++ i) + Bbox_dimensions size; + for (int i = 0; i < dimension; ++i) { - bbox_min[i] = bbox_centroid[i] - max_length; - bbox_max[i] = bbox_centroid[i] + max_length; + size[i] = (m_bbox.max)()[i] - (m_bbox.min)()[i]; } + // save orthtree attributes + m_side_per_depth.push_back(size); - Construct_point_d_from_array construct_point_d_from_array - = m_traits.construct_point_d_from_array_object(); + if constexpr (has_data) + data(root()) = m_traits.construct_root_node_contents_object()(); + } - // save orthtree attributes - m_bbox_min = construct_point_d_from_array(bbox_min); - m_side_per_depth.push_back(bbox_max[0] - bbox_min[0]); - m_root.points() = {point_range.begin(), point_range.end()}; + /*! + constructs an orthtree from a set of arguments provided to the traits constructor + */ + template = 2>> + explicit Orthtree(Args&& ... args) + : Orthtree(Traits(std::forward(args)...)) + {} + + /// copy constructor + explicit Orthtree(const Orthtree& other) : + m_traits(other.m_traits), + m_node_properties(other.m_node_properties), + m_node_contents(m_node_properties), + m_node_depths(m_node_properties.template get_property("depths")), + m_node_coordinates(m_node_properties.template get_property("coordinates")), + m_node_parents(m_node_properties.template get_property>("parents")), + m_node_children(m_node_properties.template get_property>("children")), + m_bbox(other.m_bbox), m_side_per_depth(other.m_side_per_depth) {} + + /// move constructor + explicit Orthtree(Orthtree&& other) : + m_traits(other.m_traits), + m_node_properties(std::move(other.m_node_properties)), + m_node_contents(m_node_properties), + m_node_depths(m_node_properties.template get_property("depths")), + m_node_coordinates(m_node_properties.template get_property("coordinates")), + m_node_parents(m_node_properties.template get_property>("parents")), + m_node_children(m_node_properties.template get_property>("children")), + m_bbox(other.m_bbox), m_side_per_depth(other.m_side_per_depth) + { + other.m_node_properties.emplace(); } /// @} - /// \cond SKIP_IN_MANUAL - - // copy constructor - Orthtree (const Orthtree& other) - : m_traits (other.m_traits) - , m_range (other.m_range) - , m_point_map (other.m_point_map) - , m_root (other.m_root.deep_copy()) - , m_bbox_min (other.m_bbox_min) - , m_side_per_depth(other.m_side_per_depth) - { } - - // move constructor - Orthtree (Orthtree&& other) - : m_traits (other.m_traits) - , m_range (other.m_range) - , m_point_map (other.m_point_map) - , m_root (other.m_root) - , m_bbox_min (other.m_bbox_min) - , m_side_per_depth(other.m_side_per_depth) - { - other.m_root = Node(Node(), 0); - } // Non-necessary but just to be clear on the rule of 5: - // assignment operators deleted (PointRange is a ref) - Orthtree& operator= (const Orthtree& other) = delete; - Orthtree& operator= (Orthtree&& other) = delete; - // Destructor - ~Orthtree() - { - std::queue nodes; - nodes.push(m_root); - while (!nodes.empty()) - { - Node node = nodes.front(); - nodes.pop(); - if (!node.is_leaf()) - for (std::size_t i = 0; i < Degree::value; ++ i) - nodes.push (node[i]); - node.free(); - } - } + // assignment operators deleted + Orthtree& operator=(const Orthtree& other) = delete; - // move constructor - /// \endcond + Orthtree& operator=(Orthtree&& other) = delete; /// \name Tree Building /// @{ @@ -305,105 +324,82 @@ class Orthtree /*! \brief recursively subdivides the orthtree until it meets the given criteria. - The split predicate is a `std::function` that takes a `Node` and - returns a Boolean value (where `true` implies that a `Node` needs to - be split, `false` that the `Node` should be a leaf). + The split predicate should return `true` if a leaf node should be split and `false` otherwise. This function may be called several times with different predicates: in that case, nodes already split are left unaltered, while nodes that were not split and for which `split_predicate` returns `true` are split. - \param split_predicate determines whether or not a node needs to - be subdivided. + \param split_predicate determines whether or not a leaf node needs to be subdivided. */ void refine(const Split_predicate& split_predicate) { - // If the tree has already been refined, reset it - if (!m_root.is_leaf()){ - std::queue nodes; - for (std::size_t i = 0; i < Degree::value; ++ i) - nodes.push (m_root[i]); - while (!nodes.empty()) - { - Node node = nodes.front(); - nodes.pop(); - if (!node.is_leaf()) - for (std::size_t i = 0; i < Degree::value; ++ i) - nodes.push (node[i]); - node.free(); - } - m_root.unsplit(); - } - - // Reset the side length map, too - m_side_per_depth.resize(1); - // Initialize a queue of nodes that need to be refined - std::queue todo; - todo.push(m_root); + std::queue todo; + todo.push(0); // Process items in the queue until it's consumed fully while (!todo.empty()) { // Get the next element - Node current = todo.front(); + auto current = todo.front(); todo.pop(); // Check if this node needs to be processed - if (split_predicate(current)) { + if (split_predicate(current, *this)) { - // Check if we've reached a new max depth - if (current.depth() == depth()) { + // Split the node, redistributing its contents to its children + split(current); - // Update the side length map - m_side_per_depth.push_back(*(m_side_per_depth.end() - 1) / 2); - } + } - // Split the node, redistributing its points to its children - split(current); + // Check if the node has children which need to be processed + if (!is_leaf(current)) { // Process each of its children - for (int i = 0; i < Degree::value; ++i) - todo.push(current[i]); - + for (int i = 0; i < degree; ++i) + todo.push(child(current, i)); } } } /*! - - \brief Convenience overload that refines an orthtree using a - maximum depth and maximum number of points in a node as split + \brief convenience overload that refines an orthtree using a + maximum depth and maximum number of contained elements in a node as split predicate. This is equivalent to calling - `refine(Orthtrees::Maximum_depth_and_maximum_number_of_inliers(max_depth, + `refine(Orthtrees::Maximum_depth_and_maximum_contained_elements(max_depth, bucket_size))`. The refinement is stopped as soon as one of the conditions is - violated: if a node has more inliers than `bucket_size` but is + violated: if a node contains more elements than `bucket_size` but is already at `max_depth`, it is not split. Similarly, a node that is - at a depth smaller than `max_depth` but already has fewer inliers + at a depth smaller than `max_depth` but already contains fewer elements than `bucket_size`, it is not split. + \warning This convenience method is only appropriate for trees with traits classes where + `Node_data` is a model of `Range`. `RandomAccessRange` is suggested for performance. + \param max_depth deepest a tree is allowed to be (nodes at this depth will not be split). - \param bucket_size maximum points a node is allowed to contain. + \param bucket_size maximum number of items a node is allowed to contain. */ void refine(size_t max_depth = 10, size_t bucket_size = 20) { - refine(Orthtrees::Maximum_depth_and_maximum_number_of_inliers(max_depth, bucket_size)); + refine(Orthtrees::Maximum_depth_and_maximum_contained_elements(max_depth, bucket_size)); } /*! \brief refines the orthtree such that the difference of depth between two immediate neighbor leaves is never more than 1. + + This is done only by adding nodes, nodes are never removed. */ void grade() { // Collect all the leaf nodes - std::queue leaf_nodes; - for (Node leaf : traverse(Orthtrees::Leaves_traversal())) { - // TODO: I'd like to find a better (safer) way of doing this + std::queue leaf_nodes; + for (Node_index leaf: traverse(Orthtrees::Leaves_traversal(*this))) { leaf_nodes.push(leaf); } @@ -411,42 +407,42 @@ class Orthtree while (!leaf_nodes.empty()) { // Get the next node - Node node = leaf_nodes.front(); + Node_index node = leaf_nodes.front(); leaf_nodes.pop(); // Skip this node if it isn't a leaf anymore - if (!node.is_leaf()) + if (!is_leaf(node)) continue; // Iterate over each of the neighbors for (int direction = 0; direction < 6; ++direction) { // Get the neighbor - Node neighbor = node.adjacent_node(direction); + auto neighbor = adjacent_node(node, direction); // If it doesn't exist, skip it - if (neighbor.is_null()) + if (!neighbor) continue; // Skip if this neighbor is a direct sibling (it's guaranteed to be the same depth) // TODO: This check might be redundant, if it doesn't affect performance maybe I could remove it - if (neighbor.parent() == node.parent()) + if (parent(*neighbor) == parent(node)) continue; // If it's already been split, skip it - if (!neighbor.is_leaf()) + if (!is_leaf(*neighbor)) continue; // Check if the neighbor breaks our grading rule // TODO: could the rule be parametrized? - if ((node.depth() - neighbor.depth()) > 1) { + if ((depth(node) - depth(*neighbor)) > 1) { // Split the neighbor - split(neighbor); + split(*neighbor); // Add newly created children to the queue - for (int i = 0; i < Degree::value; ++i) { - leaf_nodes.push(neighbor[i]); + for (int i = 0; i < degree; ++i) { + leaf_nodes.push(child(*neighbor, i)); } } } @@ -459,22 +455,15 @@ class Orthtree /// @{ /*! - \brief returns the root node. + * \brief provides direct read-only access to the tree traits. */ - Node root() const { return m_root; } + const Traits& traits() const { return m_traits; } /*! - \brief Convenience function to access the child nodes of the root - node by their indices. - - `my_tree[5]` is equivalent to `my_tree.root()[5]`. - - \sa `Node::operator[]()` - - \param index the index of the child node. - \return the accessed node. + \brief provides access to the root node, and by + extension the rest of the tree. */ - Node operator[](std::size_t index) const { return m_root[index]; } + Node_index root() const { return 0; } /*! \brief returns the deepest level reached by a leaf node in this tree (root being level 0). @@ -482,54 +471,156 @@ class Orthtree std::size_t depth() const { return m_side_per_depth.size() - 1; } /*! - \brief constructs a node range using a tree-traversal function. + \brief constructs a node index range using a tree-traversal function. - This method allows to iterate on the nodes of the tree with a + This method allows iteration over the nodes of the tree with a user-selected order (preorder, postorder, leaves-only, etc.). - \tparam Traversal model of `OrthtreeTraversal` that provides functions - compatible with the type of the orthree + \tparam Traversal a model of `OrthtreeTraversal` - \param traversal the instance of `Traversal` used + \param traversal class defining the traversal strategy - \return a forward input iterator over the nodes of the tree + \return a `ForwardRange` over the node indices of the tree */ - template - Node_range traverse(const Traversal &traversal = Traversal()) const { + template + Node_index_range traverse(Traversal traversal) const { - Node first = traversal.first(m_root); + Node_index first = traversal.first_index(); - Node_traversal_method_const next - = [&](const Node& n) -> Node { return traversal.next(n); }; + auto next = [=](const Self&, Node_index index) -> std::optional { + return traversal.next_index(index); + }; - return boost::make_iterator_range(Traversal_iterator(first, next), - Traversal_iterator()); + return boost::make_iterator_range(Index_traversal_iterator(*this, first, next), + Index_traversal_iterator()); } + /*! - \brief constructs the bounding box of a node. + \brief convenience method for using a traversal without constructing it yourself + + \tparam Traversal a model of `OrthtreeTraversal` - \note The object constructed is not the bounding box of the point - subset inside the node, but the bounding box of the node itself - (cubic). + \param args Arguments to to pass to the traversal's constructor, excluding the first (always an orthtree reference) + + \return a `ForwardRange` over the node indices of the tree */ - Bbox bbox(const Node &node) const { - // Determine the side length of this node - FT size = m_side_per_depth[node.depth()]; + template + Node_index_range traverse(Args&& ...args) const { + return traverse(Traversal{*this, std::forward(args)...}); + } - // Determine the location this node should be split - Array min_corner; - Array max_corner; - for (int i = 0; i < Dimension::value; i++) { + // TODO shall we document it? + FT + compute_cartesian_coordinate(std::uint32_t gc, std::size_t depth, int ci) const + { + CGAL_assertion(depth <= m_side_per_depth.size()); + // an odd coordinate will be first compute at the current depth, + // while an even coordinate has already been computed at a previous depth. + // So while the coordinate is even, we decrease the depth to end up of the first + // non-even coordinate to compute it (with particular case for bbox limits). + // Note that is depth becomes too large, we might end up with incorrect coordinates + // due to rounding errors. + if (gc == (1u << depth)) return (m_bbox.max)()[ci]; // gc == 2^node_depth + if (gc == 0) return (m_bbox.min)()[ci]; + if (gc % 2 !=0) + { + FT size = depth < m_side_per_depth.size() + ? m_side_per_depth[depth][ci] + : m_side_per_depth[depth-1][ci]/FT(2); + return (m_bbox.min)()[ci] + int(gc) * size; + } + std::size_t nd = depth; + do{ + --nd; + gc = gc >> 1; + } + while((gc&1)==0); // while even, shift + return (m_bbox.min)()[ci] + int(gc) * m_side_per_depth[nd][ci]; + } + + /*! + \brief constructs the bounding box of a node. + + \note The object constructed is not the bounding box of the node's contents, + but the bounding box of the node itself. + + \param n node to generate a bounding box for - min_corner[i] = m_bbox_min[i] + (node.global_coordinates()[i] * size); - max_corner[i] = min_corner[i] + size; + \return the bounding box of the node n + */ + Bbox bbox(Node_index n) const { + using Cartesian_coordinate = std::array; + Cartesian_coordinate min_corner, max_corner; + std::size_t node_depth = depth(n); + + for (int i = 0; i < dimension; i++) + { + min_corner[i]=compute_cartesian_coordinate(global_coordinates(n)[i], node_depth, i); + max_corner[i]=compute_cartesian_coordinate(global_coordinates(n)[i]+1, node_depth, i); } + return {std::apply(m_traits.construct_point_d_object(), min_corner), + std::apply(m_traits.construct_point_d_object(), max_corner)}; + } + + /// @} + + /// \name Custom Properties + /// @{ + + /*! + \brief gets a property for nodes, adding it if it does not already exist. + + \tparam T the type of the property to add + + \param name the name of the new property + \param default_value the default value assigned to nodes for this property + + \return pair of the property map and a Boolean which is `true` if the property needed to be created + */ + template + std::pair, bool> add_property(const std::string& name, const T default_value = T()) { + auto p = m_node_properties.get_or_add_property(name, default_value); + return std::pair, bool>(Property_map(p.first), p.second); + } - // Create the bbox - Construct_bbox_d construct_bbox - = m_traits.construct_bbox_d_object(); - return construct_bbox(min_corner, max_corner); + /*! + \brief gets a property of the nodes if it exists. + + \tparam T the type of the property to retrieve + + \param name the name of the property + + \return an optional containing the property map if it exists + */ + template + std::optional> property(const std::string& name) { + auto p = m_node_properties.template get_property_if_exists(name); + if (p) + return std::optional >(Property_map(*p)); + else + return std::nullopt; + } + + /*! + \brief returns a vector of all property names. + */ + std::vector properties() const { + return m_node_properties.properties(); + } + + /*! + \brief removes the node property from the tree. + + \tparam T the type of the property to remove + + \param property the property to be removed from the tree. + + \return true if property was a valid property of the tree. + */ + template + bool remove_property(Property_map property) { + return m_node_properties.remove_property(property.array()); } /// @} @@ -538,37 +629,39 @@ class Orthtree /// @{ /*! - \brief finds the leaf node which contains the point `p`. + \brief finds the leaf node which contains a particular point in space. - Traverses the orthtree and finds the deepest cell that has a + Traverses the orthtree and finds the leaf cell that has a domain enclosing the point passed. The point passed must be within - the region enclosed by the orthtree (bbox of the root node). + the region enclosed by the orthtree (bbox of the root node). The point is contained in the + lower cell of each direction if its coordinate is lower than the center. \param point query point. - \return the node which contains the point. + + \return the index of the node which contains the point. */ - Node locate(const Point &point) const { + Node_index locate(const Point& point) const { // Make sure the point is enclosed by the orthtree - CGAL_precondition (CGAL::do_intersect(point, bbox(m_root))); + CGAL_precondition (CGAL::do_intersect(point, bbox(root()))); // Start at the root node - auto node_for_point = m_root; + Node_index node_for_point = root(); // Descend the tree until reaching a leaf node - while (!node_for_point.is_leaf()) { + while (!is_leaf(node_for_point)) { // Find the point to split around Point center = barycenter(node_for_point); // Find the index of the correct sub-node - typename Node::Local_coordinates index; - std::size_t dimension = 0; + Local_coordinates local_coords; + std::size_t dim = 0; for (const auto& r : cartesian_range(center, point)) - index[dimension ++] = (get<0>(r) < get<1>(r)); + local_coords[dim++] = (get<0>(r) <= get<1>(r)); // Find the correct sub-node of the current node - node_for_point = node_for_point[index.to_ulong()]; + node_for_point = child(node_for_point, local_coords.to_ulong()); } // Return the result @@ -576,39 +669,91 @@ class Orthtree } /*! - \brief finds the `k` nearest neighbors of `query`. + \brief finds the `k` nearest neighbors of the point `query`. - Nearest neighbors are outputted in order of increasing distance to - `query`. + Nearest neighbors are outputted in order of increasing distance to `query`. - \tparam OutputIterator a model of `OutputIterator` that accept `Point_d` objects. - \param query a query point. - \param k the number of neighbors. - \param output the output iterator. - */ + \tparam OutputIterator a model of `OutputIterator` that accepts `GeomTraits::Node_data_element` objects. + + \param query query point + \param k number of neighbors to find + \param output output iterator + + \warning Nearest neighbor searches requires `GeomTraits` to be a model of `CollectionPartitioningOrthtreeTraits`. + */ template - OutputIterator nearest_neighbors (const Point& query, - std::size_t k, - OutputIterator output) const { - Sphere query_sphere (query, (std::numeric_limits::max)()); - return nearest_k_neighbors_in_radius(query_sphere, k, output); + auto nearest_k_neighbors(const Point& query, + std::size_t k, + OutputIterator output) const -> std::enable_if_t { + Sphere query_sphere(query, (std::numeric_limits::max)()); + CGAL_precondition(k > 0); + + return nearest_k_neighbors_within_radius(query_sphere, k, output); } /*! - \brief finds the points in `sphere`. + \brief finds the elements in the sphere `query`. - Nearest neighbors are outputted in order of increasing distance to - the center of `sphere`. + Elements are outputted in order of increasing distance to + the center of the sphere. - \tparam OutputIterator a model of `OutputIterator` that accept `Point_d` objects. - \param query query sphere. - \param output output iterator. - */ + \tparam OutputIterator a model of `OutputIterator` that accepts `GeomTraits::Node_data_element` objects. + + \param query query sphere + \param output output iterator + + \warning Nearest neighbor searches requires `GeomTraits` to be a model of `CollectionPartitioningOrthtreeTraits`. + */ template - OutputIterator nearest_neighbors (const Sphere& query, OutputIterator output) const { + auto neighbors_within_radius(const Sphere& query, OutputIterator output) const -> std::enable_if_t { + return nearest_k_neighbors_within_radius(query, (std::numeric_limits::max)(), output); + } + + /*! + \brief finds at most `k` elements within a specific radius that are + nearest to the center of the sphere `query`: if `query` does not contain + at least `k` elements, only contained elements will be returned. + + This function is useful when the user already knows how sparse the elements are, + or if they do not care about elements that are too far away. + Setting a small radius may have performance benefits. + + \tparam OutputIterator must be a model of `OutputIterator` that accepts `GeomTraits::Node_data_element` + + \param query the region to search within + \param k the number of elements to find + \param output the output iterator to add the found elements to (in order of increasing distance) + + \warning Nearest neighbor searches requires `GeomTraits` to be a model of `CollectionPartitioningOrthtreeTraits`. + */ + template + auto nearest_k_neighbors_within_radius( + const Sphere& query, + std::size_t k, + OutputIterator output + ) const -> std::enable_if_t { + CGAL_precondition(k > 0); Sphere query_sphere = query; - return nearest_k_neighbors_in_radius(query_sphere, - (std::numeric_limits::max)(), output); + + // todo: this type is over-constrained, this must be made more generic + struct Node_element_with_distance { + typename Traits::Node_data_element element; + FT distance; + }; + + // Create an empty list of elements + std::vector element_list; + if (k != (std::numeric_limits::max)()) + element_list.reserve(k); + + // Invoking the recursive function adds those elements to the vector (passed by reference) + nearest_k_neighbors_recursive(query_sphere, root(), element_list, k); + + // Add all the points found to the output + for (auto& item : element_list) + *output++ = item.element; + + return output; } /*! @@ -617,15 +762,18 @@ class Orthtree \note this function requires the function `bool CGAL::do_intersect(QueryType, Traits::Bbox_d)` to be defined. - This function finds all the intersecting nodes and returns them as const pointers. + This function finds all the intersecting leaf nodes and writes their indices to the output iterator. + + \tparam Query the primitive class (e.g., sphere, ray) + \tparam OutputIterator a model of `OutputIterator` that accepts `Node_index` types - \tparam Query the primitive class (e.g. sphere, ray) - \tparam OutputIterator a model of `OutputIterator` that accepts `Node` objects \param query the intersecting primitive. \param output output iterator. + + \return the output iterator after writing */ - template - OutputIterator intersected_nodes (const Query& query, OutputIterator output) const { + template + OutputIterator intersected_nodes(const Query& query, OutputIterator output) const { return intersected_nodes_recursive(query, root(), output); } @@ -637,14 +785,17 @@ class Orthtree /*! \brief compares the topology of the orthtree with that of `rhs`. - Trees may be considered equivalent even if they contain different points. - Equivalent trees must have the same bounding box and the same node structure. - Node structure is evaluated by comparing the root nodes using the node equality operator. + Trees may be considered equivalent even if they have different contents. + Equivalent trees must have the same root bounding box and the same node structure. + + \param rhs the other orthtree + + \return `true` if the trees have the same topology, and `false` otherwise */ - bool operator==(const Self &rhs) const { + bool operator==(const Self& rhs) const { // Identical trees should have the same bounding box - if (rhs.m_bbox_min != m_bbox_min || rhs.m_side_per_depth[0] != m_side_per_depth[0]) + if (rhs.m_bbox != m_bbox || rhs.m_side_per_depth[0] != m_side_per_depth[0]) return false; // Identical trees should have the same depth @@ -652,273 +803,613 @@ class Orthtree return false; // If all else is equal, recursively compare the trees themselves - return Node::is_topology_equal(rhs.m_root, m_root); + return is_topology_equal(*this, rhs); } /*! \brief compares the topology of the orthtree with that of `rhs`. + + \param rhs the other orthtree + + \return `false` if the trees have the same topology, and `true` otherwise */ - bool operator!=(const Self &rhs) const { + bool operator!=(const Self& rhs) const { return !operator==(rhs); } /// @} - // TODO: Document this - // TODO: Could this method name be reduced to just "center" ? - Point barycenter(const Node& node) const { + /// \name Node Access + /// @{ - // Determine the side length of this node - FT size = m_side_per_depth[node.depth()]; + /*! + \brief determines whether the node specified by index `n` is a leaf node. + */ + bool is_leaf(Node_index n) const { + return !m_node_children[n].has_value(); + } - // Determine the location this node should be split - Array bary; - std::size_t i = 0; - for (const FT& f : cartesian_range(m_bbox_min)) - { - bary[i] = FT(node.global_coordinates()[i]) * size + size / FT(2) + f; - ++ i; - } + /*! + \brief determines whether the node specified by index `n` is the root node. + */ + bool is_root(Node_index n) const { + return n == 0; + } + + /*! + \brief determines the depth of the node specified. - // Convert that location into a point - Construct_point_d_from_array construct_point_d_from_array - = m_traits.construct_point_d_from_array_object(); - return construct_point_d_from_array(bary); + The root node has depth 0, its children have depth 1, and so on. + + \param n index of the node to check. + + \return the depth of the node n within its tree. + */ + std::size_t depth(Node_index n) const { + return m_node_depths[n]; } -private: // functions : + /*! + \brief retrieves a reference to the `Node_data` associated with the node specified by `n` if + `GeomTraits` is a model of `OrthtreeTraitswithData`, and `nullptr` otherwise. + */ + std::conditional_t& data(Node_index n){ + return m_node_contents[n]; + } + + /*! + \brief retrieves a const reference to the `Node_data` associated with the node specified by `n` if + `GeomTraits` is a model of `OrthtreeTraitswithData`, and `nullptr` otherwise. + */ + std::conditional_t data(Node_index n) const{ + return m_node_contents[n]; + } - void reassign_points(Node &node, Range_iterator begin, Range_iterator end, const Point ¢er, - std::bitset coord = {}, - std::size_t dimension = 0) { + /*! + \brief retrieves the global coordinates of the node. + */ + Global_coordinates global_coordinates(Node_index n) const { + return m_node_coordinates[n]; + } + + /*! + \brief retrieves the local coordinates of the node. + */ + Local_coordinates local_coordinates(Node_index n) const { + Local_coordinates result; + for (std::size_t i = 0; i < dimension; ++i) + result[i] = global_coordinates(n)[i] & 1; + return result; + } + + /*! + \brief returns this n's parent. + + \pre `!is_root()` + + \param n index of the node to retrieve the parent of + + \return the index of the parent of node n + */ + Node_index parent(Node_index n) const { + CGAL_precondition (!is_root(n)); + return *m_node_parents[n]; + } + + /*! + \brief returns a node's `i`th child. + + \pre `!is_leaf()` + + \param n index of the node to retrieve the child of + \param i in [0, 2^D) specifying the child to retrieve + + \return the index of the `i`th child of node n + */ + Node_index child(Node_index n, std::size_t i) const { + CGAL_precondition (!is_leaf(n)); + return *m_node_children[n] + i; + } + + /*! + \brief retrieves an arbitrary descendant of the node specified by `node`. + + Convenience function to avoid the need to call `orthtree.child(orthtree.child(node, 0), 1)`. + + Each index in `indices` specifies which child to enter as descending the tree from `node` down. + Indices are evaluated in the order they appear as parameters, so + `descendant(root, 0, 1)` returns the second child of the first child of the root. + + \param node the node to descend + \param indices the integer indices specifying the descent to perform + + \return the index of the specified descendant node + */ + template + Node_index descendant(Node_index node, Indices... indices) { + return recursive_descendant(node, indices...); + } + + /*! + \brief convenience function for retrieving arbitrary nodes. + + Equivalent to `tree.descendant(tree.root(), indices...)`. + + \param indices the integer indices specifying the descent to perform, starting from the root + + \return the index of the specified node + */ + template + Node_index node(Indices... indices) { + return descendant(root(), indices...); + } + + /*! + \brief finds the next sibling in the parent of the node specified by the index `n`. + + Traverses the tree in increasing order of local index (e.g., 000, 001, 010, etc.) + + \param n the index of the node to find the sibling of + + \return the index of the next sibling of n + if n is not the last node in its parent, otherwise `std::nullopt`. + */ + const std::optional next_sibling(Node_index n) const { + + // Root node has no siblings + if (is_root(n)) return {}; + + // Find out which child this is + std::size_t local_coords = local_coordinates(n).to_ulong(); + + // The last child has no more siblings + if (int(local_coords) == degree - 1) + return {}; + + // The next sibling is the child of the parent with the following local coordinates + return child(parent(n), local_coords + 1); + } + + /*! + \brief finds the next sibling of the parent of the node specified by `n` if it exists. + + \param n the index node to find the sibling up of. + + \return The index of the next sibling of the parent of n + if n is not the root and its parent has a sibling, otherwise nothing. + */ + const std::optional next_sibling_up(Node_index n) const { + + // the root node has no next sibling up + if (n == 0) return {}; - // Root case: reached the last dimension - if (dimension == Dimension::value) { + auto up = std::optional{parent(n)}; + while (up) { - node[coord.to_ulong()].points() = {begin, end}; + if (next_sibling(*up)) return {next_sibling(*up)}; - return; + up = is_root(*up) ? std::optional{} : std::optional{parent(*up)}; } - // Split the point collection around the center point on this dimension - Range_iterator split_point = std::partition - (begin, end, - [&](const Range_type &a) -> bool { - // This should be done with cartesian iterator but it seems - // complicated to do efficiently - return (get(m_point_map, a)[int(dimension)] < center[int(dimension)]); - }); + return {}; + } + + /*! + \brief finds the leaf node reached when descending the tree and always choosing child 0. + + This is the starting point of a depth-first traversal. + + \param n the index of the node to find the deepest first child of. + + \return the index of the deepest first child of node n. + */ + Node_index deepest_first_child(Node_index n) const { + + auto first = n; + while (!is_leaf(first)) + first = child(first, 0); + + return first; + } + + /*! + \brief finds node reached when descending the tree to a depth `d` and always choosing child 0. + + Similar to `deepest_first_child()`, but does go to a fixed depth. + + \param n the index of the node to find the `d`th first child of. + \param d the depth to descend to. + + \return the index of the `d`th first child, nothing if the tree is not deep enough. + */ + std::optional first_child_at_depth(Node_index n, std::size_t d) const { + + std::queue todo; + todo.push(n); + + while (!todo.empty()) { + Node_index node = todo.front(); + todo.pop(); - // Further subdivide the first side of the split - std::bitset coord_left = coord; - coord_left[dimension] = false; - reassign_points(node, begin, split_point, center, coord_left, dimension + 1); + if (depth(node) == d) + return node; - // Further subdivide the second side of the split - std::bitset coord_right = coord; - coord_right[dimension] = true; - reassign_points(node, split_point, end, center, coord_right, dimension + 1); + if (!is_leaf(node)) + for (int i = 0; i < degree; ++i) + todo.push(child(node, i)); + } + return {}; } - void split(Node& node) { + /*! + \brief splits a node into subnodes. + + Only leaf nodes should be split. + When a node is split it is no longer a leaf node. + The full set of `degree` children are constructed automatically, and their values are set. + Contents of this node are _not_ propagated automatically, this is responsibility of the + `distribute_node_contents_object` in the traits class. + + \param n index of the node to split + */ + void split(Node_index n) { // Make sure the node hasn't already been split - CGAL_precondition (node.is_leaf()); + CGAL_precondition (is_leaf(n)); // Split the node to create children - node.split(); + using Local_coordinates = Local_coordinates; + m_node_children[n] = m_node_properties.emplace_group(degree); + for (std::size_t i = 0; i < degree; i++) { + + Node_index c = *m_node_children[n] + i; + + // Make sure the node isn't one of its own children + CGAL_assertion(n != *m_node_children[n] + i); + + Local_coordinates local_coordinates{i}; + for (int i = 0; i < dimension; i++) + m_node_coordinates[c][i] = (2 * m_node_coordinates[n][i]) + local_coordinates[i]; + m_node_depths[c] = m_node_depths[n] + 1; + m_node_parents[c] = n; + } + + // Check if we've reached a new max depth + if (depth(n) + 1 == m_side_per_depth.size()) { + // Update the side length map with the dimensions of the children + Bbox_dimensions size = m_side_per_depth.back(); + Bbox_dimensions child_size; + for (int i = 0; i < dimension; ++i) + child_size[i] = size[i] / FT(2); + m_side_per_depth.push_back(child_size); + } - // Find the point to around which the node is split - Point center = barycenter(node); + // Find the point around which the node is split + Point center = barycenter(n); - // Add the node's points to its children - reassign_points(node, node.points().begin(), node.points().end(), center); + // Add the node's contents to its children + if constexpr (has_data) + m_traits.distribute_node_contents_object()(n, *this, center); } - bool do_intersect(const Node &node, const Sphere &sphere) const { + /*! + * \brief returns the center point of a node. + * + * @param n index of the node to find the center point for + * + * @return the center point of node n + */ + Point barycenter(Node_index n) const { + std::size_t node_depth = depth(n); + // the barycenter is computed as the lower corner of the lexicographically top child node + std::array bary; + for (int i = 0; i < dimension; i++) + bary[i] = compute_cartesian_coordinate(2 * global_coordinates(n)[i]+1, node_depth+1, i); + + return std::apply(m_traits.construct_point_d_object(), bary); + } + + /*! + \brief determines whether a pair of subtrees have the same topology. + + \param lhsNode index of a node in lhsTree + \param lhsTree an Orthtree + \param rhsNode index of a node in rhsTree + \param rhsTree another Orthtree - // Create a cubic bounding box from the node - Bbox node_cube = bbox(node); + @return true if lhsNode and rhsNode have the same topology, false otherwise + */ + static bool is_topology_equal(Node_index lhsNode, const Self& lhsTree, Node_index rhsNode, const Self& rhsTree) { + + // If one node is a leaf, and the other isn't, they're not the same + if (lhsTree.is_leaf(lhsNode) != rhsTree.is_leaf(rhsNode)) + return false; - // Check for overlap between the node's box and the sphere as a box, to quickly catch some cases - // FIXME: Activating this causes slower times -// if (!do_overlap(node_cube, sphere.bbox())) -// return false; + // If both nodes are non-leaf + if (!lhsTree.is_leaf(lhsNode)) { + + // Check all the children + for (int i = 0; i < degree; ++i) { + // If any child cell is different, they're not the same + if (!is_topology_equal(lhsTree.child(lhsNode, i), lhsTree, + rhsTree.child(rhsNode, i), rhsTree)) + return false; + } + } + + return (lhsTree.global_coordinates(lhsNode) == rhsTree.global_coordinates(rhsNode)); + } + + /*! + \brief helper function for calling `is_topology_equal()` on the root nodes of two trees. + + \param lhs an Orthtree + \param rhs another Orthtree + + \return `true` if `lhs` and `rhs` have the same topology, and `false` otherwise + */ + static bool is_topology_equal(const Self& lhs, const Self& rhs) { + return is_topology_equal(lhs.root(), lhs, rhs.root(), rhs); + } + + /*! + \brief finds the directly adjacent node in a specific direction + + \pre `direction.to_ulong < 2 * dimension` + + Adjacent nodes are found according to several properties: + - adjacent nodes may be larger than the seek node, but never smaller + - a node has at most `2 * dimension` different adjacent nodes (in 3D: left, right, up, down, front, back) + - adjacent nodes are not required to be leaf nodes + + Here's a diagram demonstrating the concept for a quadtree: + + ``` + +---------------+---------------+ + | | | + | | | + | | | + | A | | + | | | + | | | + | | | + +-------+-------+---+---+-------+ + | | | | | | + | A | (S) +---A---+ | + | | | | | | + +---+---+-------+---+---+-------+ + | | | | | | + +---+---+ A | | | + | | | | | | + +---+---+-------+-------+-------+ + ``` + + - (S) : Seek node + - A : Adjacent node + + Note how the top adjacent node is larger than the seek node. The + right adjacent node is the same size, even though it contains + further subdivisions. + + This implementation returns the adjacent node if it's found. If + there is no adjacent node in that direction, it returns a null + node. + + \param n index of the node to find a neighbor of + \param direction which way to find the adjacent node relative to + this one. Each successive bit selects the direction for the + corresponding dimension: for an octree in 3D, 010 means: negative + direction in X, position direction in Y, negative direction in Z. + + \return the index of the adjacent node if it exists, nothing otherwise. + */ + std::optional adjacent_node(Node_index n, const Local_coordinates& direction) const { + + // Direction: LEFT RIGHT DOWN UP BACK FRONT + // direction: 000 001 010 011 100 101 + + // Nodes only have up to 2*dim different adjacent nodes (since boxes have 6 sides) + CGAL_precondition(direction.to_ulong() < dimension * 2); + + // The root node has no adjacent nodes! + if (is_root(n)) return {}; + + // The least significant bit indicates the sign (which side of the node) + bool sign = direction[0]; + + // The first two bits indicate the dimension/axis (x, y, z) + uint8_t dim = uint8_t((direction >> 1).to_ulong()); + + // Create an offset so that the bit-significance lines up with the dimension (e.g., 1, 2, 4 --> 001, 010, 100) + int8_t offset = (uint8_t) 1 << dim; + + // Finally, apply the sign to the offset + offset = (sign ? offset : -offset); + + // Check if this child has the opposite sign along the direction's axis + if (local_coordinates(n)[dim] != sign) { + // This means the adjacent node is a direct sibling, the offset can be applied easily! + return {child(parent(n), local_coordinates(n).to_ulong() + offset)}; + } + + // Find the parent's neighbor in that direction, if it exists + auto adjacent_node_of_parent = adjacent_node(parent(n), direction); + + // If the parent has no neighbor, then this node doesn't have one + if (!adjacent_node_of_parent) return {}; + + // If the parent's adjacent node has no children, then it's this node's adjacent node + if (is_leaf(*adjacent_node_of_parent)) + return adjacent_node_of_parent; + + // Return the nearest node of the parent by subtracting the offset instead of adding + return {child(*adjacent_node_of_parent, local_coordinates(n).to_ulong() - offset)}; + } + + /*! + \brief equivalent to `adjacent_node()`, with an adjacency direction rather than a bitset. + + \param n index of the node to find a neighbor of + \param adjacency which way to find the adjacent node relative to this one + */ + std::optional adjacent_node(Node_index n, Adjacency adjacency) const { + return adjacent_node(n, std::bitset(static_cast(adjacency))); + } + + /// @} + +private: // functions : + + Node_index recursive_descendant(Node_index node, std::size_t i) { return child(node, i); } + + template + Node_index recursive_descendant(Node_index node, std::size_t i, Indices... remaining_indices) { + return recursive_descendant(child(node, i), remaining_indices...); + } + + bool do_intersect(Node_index n, const Sphere& sphere) const { + + // Create a bounding box from the node + Bbox node_box = bbox(n); // Check for intersection between the node and the sphere - return CGAL::do_intersect(node_cube, sphere); + return CGAL::do_intersect(node_box, sphere); } - // TODO: There has to be a better way than using structs like these! - struct Point_with_distance { - Point point; - FT distance; - }; - struct Node_index_with_distance { - typename Node::Local_coordinates index; - FT distance; + template + Node_output_iterator intersected_nodes_recursive(const Query& query, Node_index node, + Node_output_iterator output) const { - Node_index_with_distance (const typename Node::Local_coordinates& index, - const FT& distance) - : index(index), distance(distance) - { } - }; + // Check if the current node intersects with the query + if (CGAL::do_intersect(query, bbox(node))) { - void nearest_k_neighbors_recursive(Sphere& search_bounds, const Node &node, - std::vector &results, FT epsilon = 0) const { + // if this node is a leaf, then it's considered an intersecting node + if (is_leaf(node)) { + *output++ = node; + return output; + } + + // Otherwise, each of the children need to be checked + for (int i = 0; i < degree; ++i) { + intersected_nodes_recursive(query, child(node, i), output); + } + } + return output; + } + + template + auto nearest_k_neighbors_recursive( + Sphere& search_bounds, + Node_index node, + std::vector& results, + std::size_t k, + FT epsilon = 0) const -> std::enable_if_t { // Check whether the node has children - if (node.is_leaf()) { + if (is_leaf(node)) { // Base case: the node has no children - // Loop through each of the points contained by the node + // Loop through each of the elements contained by the node // Note: there might be none, and that should be fine! - for (auto point_index : node.points()) { - - // Retrieve each point from the orthtree's point map - auto point = get(m_point_map, point_index); + for (auto& e : data(node)) { - // Pair that point with its distance from the search point - Point_with_distance current_point_with_distance = - {point, squared_distance(point, search_bounds.center())}; + // Pair that element with its distance from the search point + Result current_element_with_distance = + { e, m_traits.squared_distance_of_element_object()(e, m_traits.construct_center_d_object()(search_bounds)) }; - // Check if the new point is within the bounds - if (current_point_with_distance.distance < search_bounds.squared_radius()) { + // Check if the new element is within the bounds + if (current_element_with_distance.distance < m_traits.compute_squared_radius_d_object()(search_bounds)) { // Check if the results list is full - if (results.size() == results.capacity()) { - - // Delete a point if we need to make room + if (results.size() == k) { + // Delete a element if we need to make room results.pop_back(); } - // Add the new point - results.push_back(current_point_with_distance); + // Add the new element + results.push_back(current_element_with_distance); // Sort the list - std::sort(results.begin(), results.end(), [=](auto &left, auto &right) { + std::sort(results.begin(), results.end(), [=](auto& left, auto& right) { return left.distance < right.distance; - }); + }); // Check if the results list is full - if (results.size() == results.capacity()) { + if (results.size() == k) { // Set the search radius - search_bounds = Sphere(search_bounds.center(), results.back().distance + epsilon); + search_bounds = m_traits.construct_sphere_d_object()(m_traits.construct_center_d_object()(search_bounds), results.back().distance + epsilon); } } } - } else { + } + else { + + struct Node_index_with_distance { + Node_index index; + FT distance; + + Node_index_with_distance(const Node_index& index, FT distance) : + index(index), distance(distance) {} + }; // Recursive case: the node has children // Create a list to map children to their distances std::vector children_with_distances; - children_with_distances.reserve(Degree::value); + children_with_distances.reserve(Self::degree); // Fill the list with child nodes - for (int index = 0; index < Degree::value; ++index) { - Node child_node = node[index]; + for (int i = 0; i < Self::degree; ++i) { + auto child_node = child(node, i); + + FT squared_distance = 0; + Point c = m_traits.construct_center_d_object()(search_bounds); + Point b = barycenter(child_node); + for (const auto r : cartesian_range(c, b)) { + FT d = (get<0>(r) - get<1>(r)); + squared_distance += d * d; + } // Add a child to the list, with its distance - children_with_distances.emplace_back(typename Node::Local_coordinates(index), - CGAL::squared_distance(search_bounds.center(), barycenter(child_node))); + children_with_distances.emplace_back( + child_node, + squared_distance + ); } // Sort the children by their distance from the search point - std::sort(children_with_distances.begin(), children_with_distances.end(), [=](auto &left, auto &right) { + std::sort(children_with_distances.begin(), children_with_distances.end(), [=](auto& left, auto& right) { return left.distance < right.distance; - }); + }); // Loop over the children for (auto child_with_distance : children_with_distances) { - Node child_node = node[child_with_distance.index.to_ulong()]; // Check whether the bounding box of the child intersects with the search bounds - if (do_intersect(child_node, search_bounds)) { + if (CGAL::do_intersect(bbox(child_with_distance.index), search_bounds)) { // Recursively invoke this function - nearest_k_neighbors_recursive(search_bounds, child_node, results); + nearest_k_neighbors_recursive(search_bounds, child_with_distance.index, results, k); } } } } - template - Node_output_iterator intersected_nodes_recursive(const Query &query, const Node &node, - Node_output_iterator output) const { - - // Check if the current node intersects with the query - if (CGAL::do_intersect(query, bbox(node))) { - - // if this node is a leaf, than it's considered an intersecting node - if (node.is_leaf()) { - *output++ = node; - return output; - } - - // Otherwise, each of the children need to be checked - for (int i = 0; i < Degree::value; ++i) { - intersected_nodes_recursive(query, node[i], output); - } - } - return output; - } - - /*! - \brief finds the `k` points within a specific radius that are - nearest to the center of `query_sphere`. - - This function guarantees that there are no closer points than the ones returned, - but it does not guarantee that it will return at least `k` points. - For a query where the search radius encloses `k` or fewer points, all enclosed points will be returned. - If the search radius is too small, no points may be returned. - This function is useful when the user already knows how sparse the points are, - or if they do not care about points that are too far away. - Setting a small radius may have performance benefits. - - \tparam OutputIterator must be a model of `OutputIterator` that accepts points - \param query_sphere the region to search within - \param k the number of points to find - \param output the output iterator to add the found points to (in order of increasing distance) - */ - template - OutputIterator nearest_k_neighbors_in_radius - (Sphere& query_sphere, - std::size_t k, OutputIterator output) const { - - // Create an empty list of points - std::vector points_list; - if (k != (std::numeric_limits::max)()) - points_list.reserve(k); - - // Invoking the recursive function adds those points to the vector (passed by reference) - nearest_k_neighbors_recursive(query_sphere, m_root, points_list); - - // Add all the points found to the output - for (auto &item : points_list) - *output++ = item.point; - - return output; - } - public: /// \cond SKIP_IN_MANUAL - void dump_to_polylines (std::ostream& os) const - { - for (const Node& n : traverse()) - if (n.is_leaf()) - { + void dump_to_polylines(std::ostream& os) const { + for (const auto n: traverse()) + if (is_leaf(n)) { Bbox box = bbox(n); - dump_box_to_polylines (box, os); + dump_box_to_polylines(box, os); } } - void dump_box_to_polylines (const Bbox_2& box, std::ostream& os) const - { + void dump_box_to_polylines(const Bbox_2& box, std::ostream& os) const { // dump in 3D for visualisation os << "5 " << box.xmin() << " " << box.ymin() << " 0 " @@ -927,8 +1418,8 @@ class Orthtree << box.xmax() << " " << box.ymin() << " 0 " << box.xmin() << " " << box.ymin() << " 0" << std::endl; } - void dump_box_to_polylines (const Bbox_3& box, std::ostream& os) const - { + + void dump_box_to_polylines(const Bbox_3& box, std::ostream& os) const { // Back face os << "5 " << box.xmin() << " " << box.ymin() << " " << box.zmin() << " " @@ -960,17 +1451,21 @@ class Orthtree << box.xmax() << " " << box.ymax() << " " << box.zmax() << std::endl; } - friend std::ostream& operator<< (std::ostream& os, const Self& orthtree) - { - // Create a range of nodes - auto nodes = orthtree.traverse(Orthtrees::Preorder_traversal()); - // Iterate over the range - for (auto &n : nodes) { + std::string to_string(Node_index node) { + std::stringstream stream; + internal::print_orthtree_node(stream, node, *this); + return stream.str(); + } + + friend std::ostream& operator<<(std::ostream& os, const Self& orthtree) { + // Iterate over all nodes + for (auto n: orthtree.traverse(Orthtrees::Preorder_traversal(orthtree))) { // Show the depth - for (int i = 0; i < n.depth(); ++i) + for (std::size_t i = 0; i < orthtree.depth(n); ++i) os << ". "; // Print the node - os << n << std::endl; + internal::print_orthtree_node(os, n, orthtree); + os << std::endl; } return os; } @@ -981,6 +1476,4 @@ class Orthtree } // namespace CGAL -#include - #endif // CGAL_ORTHTREE_H diff --git a/Orthtree/include/CGAL/Orthtree/Cartesian_ranges.h b/Orthtree/include/CGAL/Orthtree/Cartesian_ranges.h index 38b4ef3f5aa0..80cdc3cd3e5e 100644 --- a/Orthtree/include/CGAL/Orthtree/Cartesian_ranges.h +++ b/Orthtree/include/CGAL/Orthtree/Cartesian_ranges.h @@ -29,8 +29,8 @@ namespace internal template struct Cartesian_ranges { - typedef typename Traits::Point_d Point; - typedef typename Traits::Cartesian_const_iterator_d Cartesian_const_iterator; + using Point = typename Traits::Point_d; + using Cartesian_const_iterator = typename Traits::Cartesian_const_iterator_d; using Range_single = CGAL::Iterator_range; diff --git a/Orthtree/include/CGAL/Orthtree/IO.h b/Orthtree/include/CGAL/Orthtree/IO.h index 3c5748397c9d..1bf6bdd71438 100644 --- a/Orthtree/include/CGAL/Orthtree/IO.h +++ b/Orthtree/include/CGAL/Orthtree/IO.h @@ -22,46 +22,39 @@ namespace CGAL namespace internal { -template -std::ostream& print_orthtree_node(std::ostream& os, const Node& node) +template +std::ostream& print_orthtree_node(std::ostream& os, const Node_index& node, const Tree &tree) { - // Show the depth of the node -// for (int i = 0; i < node.depth(); ++i) -// os << ". "; // Wrap information in brackets os << "{ "; // Index identifies which child this is os << "("; - for (std::size_t i = 0; i < node.local_coordinates().size(); ++ i) - os << node.local_coordinates()[i]; + for (std::size_t i = 0; i < tree.local_coordinates(node).size(); ++ i) + os << tree.local_coordinates(node)[i]; os << ") "; // Location os << "( "; - for (const auto& b : node.global_coordinates()) + for (const auto& b : tree.global_coordinates(node)) os << b << " "; os << ") "; // Depth os << "(" - << +node.depth() // The + forces printing as an int instead of a char + << +tree.depth(node) // The + forces printing as an int instead of a char << ") "; os << "(" - << node.size() + << tree.data(node).size() << ") "; -// // If a node has points, indicate how many -// if (!node.is_empty()) -// os << "[" << node.num_points() << " points] "; - // If a node is a leaf, mark it - os << (node.is_leaf() ? "[leaf] " : ""); + os << (tree.is_leaf(node) ? "[leaf] " : ""); // If a node is root, mark it - os << (node.is_root() ? "[root] " : ""); + os << (tree.is_root(node) ? "[root] " : ""); // Wrap information in brackets os << "}"; diff --git a/Orthtree/include/CGAL/Orthtree/Node.h b/Orthtree/include/CGAL/Orthtree/Node.h deleted file mode 100644 index 777066a295d8..000000000000 --- a/Orthtree/include/CGAL/Orthtree/Node.h +++ /dev/null @@ -1,619 +0,0 @@ -// Copyright (c) 2007-2020 INRIA (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) : Jackson Campolattaro, Cédric Portaneri, Tong Zhao - -#ifndef CGAL_ORTHTREE_NODE_H -#define CGAL_ORTHTREE_NODE_H - -#include - -#include - -#include - -#include -#include -#include -#include -#include - -namespace CGAL { - -/// \cond SKIP_IN_MANUAL -namespace Orthtrees -{ - -// Non-documented, or testing purpose only -struct Node_access -{ - template - static Node create_node (Node parent, LC local_coordinates) - { - return Node(parent, local_coordinates); - } - - template - static typename Node::Point_range& points(Node node) { return node.points(); } - - template - static void split(Node node) { return node.split(); } - - template - static void free(Node node) - { - typedef Dimension_tag<(2 << (Node::Dimension::value - 1))> Degree; - std::queue nodes; - nodes.push(node); - while (!nodes.empty()) - { - Node node = nodes.front(); - nodes.pop(); - if (!node.is_leaf()){ - for (std::size_t i = 0; i < Degree::value; ++ i){ - nodes.push (node[i]); - } - } - node.free(); - } - } - -}; - -} // namespace Orthtrees -/// \endcond - -/*! - - \brief represents a single node of the tree. Alternatively referred - to as a cell, orthant, or sub-tree. - - A `Node` is a lightweight object and thus generally passed by - copy. It is also a model of `ConstRange` with value type `Traits::Point_d`. - - \cgalModels{ConstRange} - */ -template -class Orthtree::Node -{ - -public: - - /// \name Types - /// @{ - - typedef Orthtree Enclosing; ///< Orthtree type (enclosing class). - typedef typename Enclosing::Dimension Dimension; ///< Dimension type. - typedef typename Enclosing::Degree Degree; ///< Degree type. - - /*! - \brief Self typedef for convenience. - */ - typedef typename Orthtree::Node Self; - - - /// \cond SKIP_IN_MANUAL - /*! - * \brief Array for containing the child nodes of this node. - */ - typedef std::array Children; - /// \endcond - - /*! - \brief Set of bits representing this node's relationship to its parent. - - Equivalent to an array of Booleans, where index[0] is whether `x` - is greater, index[1] is whether `y` is greater, index[2] is whether - `z` is greater, and so on for higher dimensions if needed. - Used to represent a node's relationship to the center of its parent. - */ - typedef std::bitset Local_coordinates; - - /*! - \brief Coordinates representing this node's relationship - with the rest of the tree. - - Each value `(x, y, z, ...)` of global coordinates is calculated by doubling - the parent's global coordinates and adding the local coordinates. - */ - typedef std::array Global_coordinates; - - - typedef typename PointRange::const_iterator const_iterator; ///< constant iterator type. - - /*! - \brief Adjacency directions. - */ - typedef typename Traits::Adjacency Adjacency; - - /// @} - -private: - - /// \cond SKIP_IN_MANUAL - /*! - \brief point range type. - */ - typedef typename PointRange::iterator iterator; ///< constant iterator type. - typedef boost::iterator_range Point_range; - /// \endcond - - // make Node trivially copiabled - struct Data - { - Point_range points; - Self parent; - std::uint8_t depth; - Global_coordinates global_coordinates; - std::unique_ptr children; - - Data (Self parent) - : parent (parent), depth (0) { } - }; - - Data* m_data; - - - /// \cond SKIP_IN_MANUAL - - // Only the Orthtree class has access to the non-default - // constructor, mutators, etc. - friend Enclosing; - - // Hidden class to access methods for testing purposes - friend Orthtrees::Node_access; - - /*! - * \brief Access to the content held by this node - * \return a reference to the collection of point indices - */ - Point_range &points() { return m_data->points; } - const Point_range &points() const { return m_data->points; } - - /// \name Construction - /// @{ - - /*! - \brief creates a new node, optionally as the child of a parent - - If no parent is provided, the node created is assumed to be the - root of a tree. This means that `parent.is_null()` returns - `true`, and the depth is zero. If a parent is provided, the node - becomes the child of that parent. In that case, an index should - be passed, telling this node its relationship to its parent. - Depth and global coordinates are automatically determined in the - constructor, and should generally be considered immutable after - construction. - - \param parent the node containing this one - \param index this node's relationship to its parent - */ - explicit Node(Self parent, Local_coordinates local_coordinates) - : m_data (new Data(parent)) { - - if (!parent.is_null()) { - - m_data->depth = parent.m_data->depth + 1; - - for (int i = 0; i < Dimension::value; i++) - m_data->global_coordinates[i] = (2 * parent.m_data->global_coordinates[i]) + local_coordinates[i]; - - } - else - for (int i = 0; i < Dimension::value; i++) - m_data->global_coordinates[i] = 0; - } - - void free() { delete m_data; } - - Node deep_copy(Self parent = Node()) const - { - if (is_null()) - return Node(); - - Node out; - out.m_data = new Data(parent); - - out.m_data->points = m_data->points; - out.m_data->depth = m_data->depth; - out.m_data->global_coordinates = m_data->global_coordinates; - std::unique_ptr children; - if (!is_leaf()) - { - out.m_data->children = std::make_unique(); - for (int index = 0; index < Degree::value; index++) - (*out.m_data->children)[index] = (*this)[index].deep_copy(out); - } - return out; - } - - /// @} - - /// \name Mutators - /// @{ - - /*! - \brief splits a node into subnodes. - - Only leaf nodes should be split. - When a node is split it is no longer a leaf node. - A number of `Degree::value` children are constructed automatically, and their values are set. - Contents of this node are _not_ propagated automatically. - It is the responsibility of the caller to redistribute the points contained by a node after splitting - */ - void split() { - - CGAL_precondition (is_leaf()); - - m_data->children = std::make_unique(); - for (int index = 0; index < Degree::value; index++) { - - (*m_data->children)[index] = std::move(Self(*this, {Local_coordinates(index)})); - - } - } - - /*! - * \brief eliminates this node's children, making it a leaf node. - * - * When a node is un-split, its children are automatically deleted. - * After un-splitting a node it will be considered a leaf node. - */ - void unsplit() { - - m_data->children.reset(); - } - - /// @} - - /// \endcond - - -public: - - /// \cond SKIP_IN_MANUAL - // Default creates null node - Node() : m_data(nullptr) { } - - // Comparison operator - bool operator< (const Node& other) const { return m_data < other.m_data; } - /// \endcond - - /// \name Type & Location - /// @{ - - /*! - \brief returns `true` if the node is null, `false` otherwise. - */ - bool is_null() const { return (m_data == nullptr); } - - /*! - \brief returns `true` if the node has no parent, `false` otherwise. - \pre `!is_null()` - */ - bool is_root() const - { - CGAL_precondition(!is_null()); - return m_data->parent.is_null(); - } - - /*! - \brief returns `true` if the node has no children, `false` otherwise. - \pre `!is_null()` - */ - bool is_leaf() const - { - CGAL_precondition(!is_null()); - return (!m_data->children); - } - - /*! - \brief returns this node's depth. - \pre `!is_null()` - */ - std::uint8_t depth() const - { - CGAL_precondition (!is_null()); - return m_data->depth; - } - - /*! - \brief returns this node's local coordinates (in relation to its parent). - \pre `!is_null()` - */ - Local_coordinates local_coordinates() const { - - CGAL_precondition (!is_null()); - // TODO: There must be a better way of doing this! - - Local_coordinates result; - - for (std::size_t i = 0; i < Dimension::value; ++ i) - result[i] = global_coordinates()[i] & 1; - - return result; - } - - /*! - \brief returns this node's global coordinates. - \pre `!is_null()` - */ - Global_coordinates global_coordinates() const - { - CGAL_precondition (!is_null()); - return m_data->global_coordinates; - } - - - /// \name Adjacency - /// @{ - - /*! - \brief returns this node's parent. - \pre `!is_null()` - */ - Self parent() const - { - CGAL_precondition (!is_null()); - return m_data->parent; - } - - /*! - \brief returns the nth child of this node. - - \pre `!is_null()` - \pre `!is_leaf()` - \pre `0 <= index && index < Degree::value` - - The orthtree subdivides the space in 2 on each dimension - available, so a child node can be accessed by selecting a Boolean - value for each dimension. The `index` parameter is thus - interpreted as a bitmap, where each bit matches a dimension - (starting by the least significant bit for coordinate X). - - For example, in the case of an octree (dimension 3): - - - index 0 (000 in binary) is the children on the "minimum corner" (xmin, ymin, zmin) - - index 1 (001 in binary) is on (xmax, ymin, zmin) - - index 2 (010 in binary) is on (xmin, ymax, zmin) - - index 6 (101 in binary) is on (xmax, ymin, zmax) - - Diagram of a quadtree: - - ``` - Children of the current node: - - Y axis - A - | +-------------------+-------------------+ - | | Coord: Ymax Xmin | Coord: Ymax Xmax | - | | Bitmap: 1 0 | Bitmap: 1 1 | - | | | | - | | -> index = 2 | -> index = 3 | - | | | | - | | | | - | | | | - | | | | - | +-------------------+-------------------+ - | | Coord: Ymin Xmin | Coord: Ymin Xmax | - | | Bitmap: 0 0 | Bitmap: 0 1 | - | | | | - | | -> index = 0 | -> index = 1 | - | | | | - | | | | - | | | | - | | | | - | +-------------------+-------------------+ - | - +--------------------------------------------> X axis - 0 - ``` - - The operator can be chained. For example, `n[5][2][3]` returns the - third child of the second child of the fifth child of a node `n`. - */ - Self operator[](std::size_t index) const { - - CGAL_precondition (!is_null()); - CGAL_precondition (!is_leaf()); - CGAL_precondition (index < Degree::value); - - return (*m_data->children)[index]; - } - - /*! - \brief finds the directly adjacent node in a specific direction - - \pre `!is_null()` - \pre `direction.to_ulong < 2 * Dimension::value` - - Adjacent nodes are found according to several properties: - - adjacent nodes may be larger than the seek node, but never smaller - - a node has at most `2 * Dimension::value` different adjacent nodes (in 3D: left, right, up, down, front, back) - - adjacent nodes are not required to be leaf nodes - - Here's a diagram demonstrating the concept for a Quadtree: - - ``` - +---------------+---------------+ - | | | - | | | - | | | - | A | | - | | | - | | | - | | | - +-------+-------+---+---+-------+ - | | | | | | - | A | (S) +---A---+ | - | | | | | | - +---+---+-------+---+---+-------+ - | | | | | | - +---+---+ A | | | - | | | | | | - +---+---+-------+-------+-------+ - ``` - - - (S) : Seek node - - A : Adjacent node - - Note how the top adjacent node is larger than the seek node. The - right adjacent node is the same size, even though it contains - further subdivisions. - - This implementation returns the adjacent node if it's found. If - there is no adjacent node in that direction, it returns a null - node. - - \param direction which way to find the adjacent node relative to - this one. Each successive bit selects the direction for the - corresponding dimension: for an Octree in 3D, 010 means: negative - direction in X, position direction in Y, negative direction in Z. - - \return the adjacent node if it exists, a null node otherwise. - */ - Self adjacent_node (Local_coordinates direction) const - { - CGAL_precondition(!is_null()); - - // Direction: LEFT RIGHT DOWN UP BACK FRONT - // direction: 000 001 010 011 100 101 - - // Nodes only have up to 2*dim different adjacent nodes (since cubes have 6 sides) - CGAL_precondition(direction.to_ulong() < Dimension::value * 2); - - // The root node has no adjacent nodes! - if (is_root()) - return Self(); - - // The least significant bit indicates the sign (which side of the node) - bool sign = direction[0]; - - // The first two bits indicate the dimension/axis (x, y, z) - uint8_t dimension = uint8_t((direction >> 1).to_ulong()); - - // Create an offset so that the bit-significance lines up with the dimension (e.g. 1, 2, 4 --> 001, 010, 100) - int8_t offset = (uint8_t) 1 << dimension; - - // Finally, apply the sign to the offset - offset = (sign ? offset : -offset); - - // Check if this child has the opposite sign along the direction's axis - if (local_coordinates()[dimension] != sign) { - - // This means the adjacent node is a direct sibling, the offset can be applied easily! - return parent()[local_coordinates().to_ulong() + offset]; - } - - // Find the parent's neighbor in that direction if it exists - Self adjacent_node_of_parent = parent().adjacent_node(direction); - - // If the parent has no neighbor, then this node doesn't have one - if (adjacent_node_of_parent.is_null()) - return Node(); - - // If the parent's adjacent node has no children, then it's this node's adjacent node - if (adjacent_node_of_parent.is_leaf()) - return adjacent_node_of_parent; - - // Return the nearest node of the parent by subtracting the offset instead of adding - return adjacent_node_of_parent[local_coordinates().to_ulong() - offset]; - - } - - /*! - \brief equivalent to `adjacent_node()`, with an adjacency direction - rather than a bitset. - */ - Self adjacent_node(Adjacency adjacency) const { - return adjacent_node(std::bitset(static_cast(adjacency))); - } - - /// @} - - /// \name Point Range - /// @{ - - /*! - \brief checks whether the node is empty of points or not. - */ - bool empty() const { - return m_data->points.empty(); - } - - /*! - \brief returns the number of points of this node. - */ - std::size_t size() const { - return std::size_t(std::distance(m_data->points.begin(), m_data->points.end())); - } - - /*! - \brief returns the iterator at the start of the collection of - points held by this node. - */ - const_iterator begin() const { return m_data->points.begin(); } - - /*! - \brief returns the iterator at the end of the collection of - points held by this node. - */ - const_iterator end() const { return m_data->points.end(); } - - /// @} - - /// \cond SKIP_IN_MANUAL - /// \name Operators - /// @{ - - /*! - * \brief compares the topology of this node to another node. - * - * \todo - * - * \param rhs node to compare with - * \return whether the nodes have different topology. - */ - bool operator==(const Self &rhs) const { - return m_data == rhs.m_data; - } - - static bool is_topology_equal (const Self& a, const Self& b) - { - CGAL_assertion (!a.is_null() && !b.is_null()); - - // If one node is a leaf, and the other isn't, they're not the same - if (a.is_leaf() != b.is_leaf()) - return false; - - // If both nodes are non-leaf nodes - if (!a.is_leaf()) { - - // Check all the children - for (int i = 0; i < Degree::value; ++i) { - // If any child cell is different, they're not the same - if (!is_topology_equal(a[i], b[i])) - return false; - } - } - - // If both nodes are leaf nodes, they must be in the same location - return (a.global_coordinates() == b.global_coordinates()); - } - - friend std::ostream& operator<< (std::ostream& os, const Self& node) - { - return internal::print_orthtree_node(os, node); - } - /// \endcond -}; - -} - -#endif //CGAL_ORTHTREE_NODE_H diff --git a/Orthtree/include/CGAL/Orthtree/Split_predicates.h b/Orthtree/include/CGAL/Orthtree/Split_predicates.h index e8d93415d5e5..c3b0d929b62e 100644 --- a/Orthtree/include/CGAL/Orthtree/Split_predicates.h +++ b/Orthtree/include/CGAL/Orthtree/Split_predicates.h @@ -18,33 +18,40 @@ namespace CGAL { +template +class Orthtree; + namespace Orthtrees { /*! \ingroup PkgOrthtreeSplitPredicates - \brief A class used to choose when a node should be split depending on the number of inliers. + \brief A class used to choose when a node should be split depending on the number of contained elements. This is a bucket size predicate that considers a node should be split if it contains more than a certain number of items. + + \warning This split predicate is only appropriate for trees with traits classes which are models of `OrthtreeTraitsWithData` + and where `Node_data` is a model of `Range`. `RandomAccessRange` is suggested for performance. */ -class Maximum_number_of_inliers { +class Maximum_contained_elements { std::size_t m_bucket_size; public: /*! - \brief creates a predicate based on the number of inliers (bucket size). + \brief creates a predicate based on the number of contained elements. */ - Maximum_number_of_inliers(std::size_t bucket_size) : + Maximum_contained_elements(std::size_t bucket_size) : m_bucket_size(bucket_size) {} /*! - \brief returns `true` if `n` should be split, `false` otherwise. + \brief returns `true` if the node with index `i` should be split, `false` otherwise. */ - template - bool operator()(const Node &n) const { - return (n.size() > m_bucket_size); + template + bool operator()(typename Orthtree::Node_index i, const Orthtree &tree) const { + return (tree.data(i).size() > m_bucket_size); } + }; /*! @@ -65,30 +72,33 @@ class Maximum_depth { Maximum_depth(std::size_t max_depth) : m_max_depth(max_depth) {} /*! - \brief returns `true` if `n` should be split, `false` otherwise. + \brief returns `true` if the node with index `i` should be split, `false` otherwise. */ - template - bool operator()(const Node &n) const { - return n.depth() < m_max_depth; + template + bool operator()(typename Orthtree::Node_index i, const Orthtree &tree) const { + return (tree.depth(i) < m_max_depth); } + }; /*! \ingroup PkgOrthtreeSplitPredicates - \brief A class used to choose when a node should be split depending on the depth and the number of inliers. + \brief A class used to choose when a node should be split depending on the depth and the number of contained elements. This predicate makes a note split if it contains more than a certain number of items and if its depth is lower than a certain limit. The refinement is stopped as soon as one of the conditions is - violated: if a node has more inliers than `bucket_size` but is + violated: if a node has more elements than `bucket_size` but is already at `max_depth`, it is not split. Similarly, a node that is - at a depth smaller than `max_depth` but already has fewer inliers + at a depth smaller than `max_depth` but already has fewer elements than `bucket_size`, it is not split. + \warning This split predicate is only appropriate for trees with traits classes which are models of `OrthtreeTraitsWithData` + and where `Node_data` is a model of `Range`. `RandomAccessRange` is suggested for performance. */ -class Maximum_depth_and_maximum_number_of_inliers { +class Maximum_depth_and_maximum_contained_elements { std::size_t m_max_depth, m_bucket_size; @@ -96,18 +106,19 @@ class Maximum_depth_and_maximum_number_of_inliers { /*! \brief creates a predicate using maximum depth or bucket size. */ - Maximum_depth_and_maximum_number_of_inliers(std::size_t max_depth, std::size_t bucket_size) : + Maximum_depth_and_maximum_contained_elements(std::size_t max_depth, std::size_t bucket_size) : m_max_depth(max_depth), m_bucket_size(bucket_size) {} /*! - \brief returns `true` if `n` should be split, `false` otherwise. + \brief returns `true` if the node with index `i` should be split, `false` otherwise. */ - template - bool operator()(const Node &n) const { - std::size_t num_points = n.size(); - std::size_t depth = n.depth(); + template + bool operator()(typename Orthtree::Node_index i, const Orthtree &tree) const { + std::size_t num_points = tree.data(i).size(); + std::size_t depth = tree.depth(i); return (num_points > m_bucket_size && depth < m_max_depth); } + }; } // Orthtrees diff --git a/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h b/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h index d8a594e4ef03..101522875380 100644 --- a/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h +++ b/Orthtree/include/CGAL/Orthtree/Traversal_iterator.h @@ -14,7 +14,10 @@ #include +#include + #include +#include #include /// \cond SKIP_IN_MANUAL @@ -22,18 +25,21 @@ namespace CGAL { /*! - * \ingroup PkgOrthtreeClasses + * \ingroup PkgOrthtreeTraversal * - * \brief + * \brief Wraps a traversal definition to produce an iterator which traverses the tree when incremented. * * \todo * - * \tparam Value + * \tparam Tree The orthtree type to iterate over */ -template -class Traversal_iterator : - public boost::iterator_facade, Value, boost::forward_traversal_tag> { - +template +class Index_traversal_iterator : public boost::iterator_facade< + Index_traversal_iterator, + const typename Tree::Node_index, + boost::forward_traversal_tag, + const typename Tree::Node_index +> { public: /// \name Types @@ -44,7 +50,9 @@ class Traversal_iterator : * * \todo */ - typedef std::function Traversal_function; + using Traversal_function = std::function(const Tree&, std::size_t)>; + + using Node_index = typename Tree::Node_index; /// @} @@ -54,44 +62,51 @@ class Traversal_iterator : /// @{ /*! - * \brief + * \brief Default constructor, creates an end sentinel * * \todo */ - Traversal_iterator() : m_value(), m_next() {} + Index_traversal_iterator() : m_next() {} /*! * \brief * * \todo * + * \param tree * \param first * \param next */ - Traversal_iterator(Value first, const Traversal_function &next) : m_value(first), m_next(next) {} + Index_traversal_iterator(const Tree& tree, Node_index first, const Traversal_function& next) : + m_tree(&tree), m_index(first), m_next(next) {} /// @} private: + friend class boost::iterator_core_access; - bool equal(Traversal_iterator const &other) const { - return m_value == other.m_value; + bool equal(Index_traversal_iterator const& other) const { + return m_index == other.m_index; } void increment() { - m_value = m_next(m_value); + // invoking increment on the sentinel is undefined behavior + m_index = m_next(*m_tree, *m_index); } - Value &dereference() const { - return const_cast(m_value); + Node_index dereference() const { + return *m_index; } private: - Value m_value; + const Tree* m_tree = nullptr; + std::optional m_index; Traversal_function m_next; + }; + } /// \endcond diff --git a/Orthtree/include/CGAL/Orthtree/Traversals.h b/Orthtree/include/CGAL/Orthtree/Traversals.h index 1ab9b77cceae..acd1db941f61 100644 --- a/Orthtree/include/CGAL/Orthtree/Traversals.h +++ b/Orthtree/include/CGAL/Orthtree/Traversals.h @@ -22,159 +22,80 @@ namespace CGAL { -/// \cond SKIP_IN_MANUAL -// Forward declaration -template -class Orthtree; -/// \endcond - namespace Orthtrees { -/// \cond SKIP_IN_MANUAL - -template -Node next_sibling(Node n) { - - // Passing null returns the first node - if (n.is_null()) - return Node(); - - // If this node has no parent, it has no siblings - if (n.parent().is_null()) - return Node(); - - // Find out which child this is - std::size_t index = n.local_coordinates().to_ulong(); - - constexpr static int degree = Node::Degree::value; - // Return null if this is the last child - if (int(index) == degree - 1) - return Node(); - - // Otherwise, return the next child - return n.parent()[index + 1]; -} - -template -Node next_sibling_up(Node n) { - - if (n.is_null()) - return Node(); - - Node up = n.parent(); - - while (!up.is_null()) { - - if (!next_sibling(up).is_null()) - return next_sibling(up); - - up = up.parent(); - } - - return Node(); -} - -template -Node deepest_first_child(Node n) { - - if (n.is_null()) - return Node(); - - // Find the deepest child on the left - Node first = n; - while (!first.is_leaf()) - first = first[0]; - return first; -} - -template -Node first_child_at_depth(Node n, std::size_t depth) { - - if (n.is_null()) - return Node(); - - std::stack todo; - todo.push(n); - - if (n.depth() == depth) - return n; - - while (!todo.empty()) - { - Node node = todo.top(); - todo.pop(); - - if (node.depth() == depth) - return node; - - if (!node.is_leaf()) - for (int i = 0; i < Node::Degree::value; ++ i) - todo.push(node[std::size_t(Node::Degree::value - 1 - i)]); - } - - return Node(); -} - -/// \endcond - /*! \ingroup PkgOrthtreeTraversal \brief A class used for performing a preorder traversal. + \tparam Tree an instance of `Orthtree` + A preorder traversal starts from the root towards the leaves. \cgalModels{OrthtreeTraversal} */ +template struct Preorder_traversal { +private: + + const Tree& m_orthtree; + +public: - template - Node first(Node root) const { - return root; + using Node_index = typename Tree::Node_index; + + Preorder_traversal(const Tree& orthtree) : m_orthtree(orthtree) {} + + Node_index first_index() const { + return m_orthtree.root(); } - template - Node next(Node n) const { + std::optional next_index(Node_index n) const { - if (n.is_leaf()) { + if (m_orthtree.is_leaf(n)) { - Node next = next_sibling(n); + auto next = m_orthtree.next_sibling(n); - if (next.is_null()) - return next_sibling_up(n); + if (!next) + next = m_orthtree.next_sibling_up(n); return next; + } else { + return m_orthtree.child(n, 0); } - else // Return the first child of this node - return n[0]; } + }; /*! \ingroup PkgOrthtreeTraversal \brief A class used for performing a postorder traversal. + \tparam Tree an instance of `Orthtree` + A postorder traversal starts from the leaves towards the root. \cgalModels{OrthtreeTraversal} */ +template struct Postorder_traversal { +private: - template - Node first(Node root) const { + const Tree& m_orthtree; - return deepest_first_child(root); - } +public: - template - Node next(Node n) const { + using Node_index = typename Tree::Node_index; - Node next = deepest_first_child(next_sibling(n)); + Postorder_traversal(const Tree& orthtree) : m_orthtree(orthtree) {} - if (next.is_null()) - next = n.parent(); + Node_index first_index() const { + return m_orthtree.deepest_first_child(m_orthtree.root()); + } - return next; + std::optional next_index(Node_index n) const { + return m_orthtree.index(next(&m_orthtree[n])); } }; @@ -182,27 +103,37 @@ struct Postorder_traversal { \ingroup PkgOrthtreeTraversal \brief A class used for performing a traversal on leaves only. - All non-leave nodes are ignored. + \tparam Tree an instance of `Orthtree` + + All non-leaf nodes are ignored. \cgalModels{OrthtreeTraversal} */ +template struct Leaves_traversal { +private: - template - Node first(Node root) const { + const Tree& m_orthtree; - return deepest_first_child(root); +public: + + using Node_index = typename Tree::Node_index; + + Leaves_traversal(const Tree& orthtree) : m_orthtree(orthtree) {} + + Node_index first_index() const { + return m_orthtree.deepest_first_child(m_orthtree.root()); } - template - Node next(Node n) const { + std::optional next_index(Node_index n) const { - Node next = deepest_first_child(next_sibling(n)); + if (m_orthtree.next_sibling(n)) + return m_orthtree.deepest_first_child(*m_orthtree.next_sibling(n)); - if (next.is_null()) - next = deepest_first_child(next_sibling_up(n)); + if (m_orthtree.next_sibling_up(n)) + return m_orthtree.deepest_first_child(*m_orthtree.next_sibling_up(n)); - return next; + return {}; } }; @@ -210,52 +141,58 @@ struct Leaves_traversal { \ingroup PkgOrthtreeTraversal \brief A class used for performing a traversal of a specific depth level. - All trees at another depth are ignored. If the selected depth is + \tparam Tree an instance of `Orthtree` + + All tree nodes at another depth are ignored. If the selected depth is + All tree nodes at another depth are ignored. If the selected depth is higher than the maximum depth of the orthtree, no node will be traversed. \cgalModels{OrthtreeTraversal} */ +template struct Level_traversal { - private: - const std::size_t depth; + const Tree& m_orthtree; + const std::size_t m_depth; public: + using Node_index = typename Tree::Node_index; + /*! constructs a `depth`-level traversal. */ - Level_traversal (std::size_t depth) : depth(depth) { } + Level_traversal(const Tree& orthtree, std::size_t depth) : m_orthtree(orthtree), m_depth(depth) {} - template - Node first(Node root) const { - return first_child_at_depth(root, depth); + Node_index first_index() const { + // assumes the tree has at least one child at m_depth + return *m_orthtree.first_child_at_depth(m_orthtree.root(), m_depth); } - template - Node next(Node n) const { - std::cerr << depth << " "; - Node next = next_sibling(n); - - if (next.is_null()) - { - Node up = n; - do - { - up = next_sibling_up(up); - if (up.is_null()) - return Node(); - next = first_child_at_depth(up, depth); - } - while (next.is_null()); + std::optional next_index(Node_index n) const { + + auto next = m_orthtree.next_sibling(n); + + if (!next) { + + auto up = n; + do { + + if (!m_orthtree.next_sibling_up(up)) + return {}; + + up = *m_orthtree.next_sibling_up(up); + next = m_orthtree.first_child_at_depth(up, m_depth); + + } while (!next); } return next; } }; -} // Orthtree +} // Orthtrees } // CGAL #endif //CGAL_ORTHTREE_TRAVERSALS_H diff --git a/Orthtree/include/CGAL/Orthtree_traits.h b/Orthtree/include/CGAL/Orthtree_traits.h new file mode 100644 index 000000000000..cc1665622e7b --- /dev/null +++ b/Orthtree/include/CGAL/Orthtree_traits.h @@ -0,0 +1,61 @@ +// Copyright (c) 2023 INRIA (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) : Sven Oesau + +#ifndef ORTHTREE_TESTS_ORTHTREE_TRAITS_H +#define ORTHTREE_TESTS_ORTHTREE_TRAITS_H + +#include + +#include +#include +#include +#include + +#include + +namespace CGAL { + +/*! + \ingroup PkgOrthtreeTraits + + Traits class for defining an orthtree using the class `CGAL::Orthtree` without storing data in the nodes. + + \tparam GeomTraits model of `Kernel`. + \tparam dimension the dimension of the ambient Euclidean space. + + \cgalModels{OrthtreeTraits} + \sa `CGAL::Octree` + \sa `CGAL::Quadtree` + \sa `CGAL::Orthtree_traits_base` +*/ +template +struct Orthtree_traits : public Orthtree_traits_base { +public: + using Base = Orthtree_traits_base; + using Self = Orthtree_traits; + using Tree = Orthtree; + + using Node_index = typename Base::Node_index; + + Orthtree_traits() {} + + auto construct_root_node_bbox_object() const { + return [&]() -> typename Self::Bbox_d { + return {std::apply(Self::construct_point_d_object(), std::array{-1.0, -1.0, -1.0}), + std::apply(Self::construct_point_d_object(), std::array{1.0, 1.0, 1.0})}; + }; + } +}; + +} + + +#endif //ORTHTREE_TESTS_ORTHTREE_TRAITS_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_2.h b/Orthtree/include/CGAL/Orthtree_traits_2.h deleted file mode 100644 index b849ae450f5f..000000000000 --- a/Orthtree/include/CGAL/Orthtree_traits_2.h +++ /dev/null @@ -1,143 +0,0 @@ -// Copyright (c) 2020 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) : Simon Giraudot - -#ifndef CGAL_ORTHTREE_TRAITS_2_H -#define CGAL_ORTHTREE_TRAITS_2_H - -#include - -#include -#include - -namespace CGAL -{ - -/*! - \ingroup PkgOrthtreeTraits - - The class `Orthtree_traits_2` can be used as a template parameter of - the `Orthtree` class. - - \tparam GeomTraits model of `Kernel`. - - \cgalModels{OrthtreeTraits} - \sa `CGAL::Quadtree` - \sa `CGAL::Orthtree_traits_3` - \sa `CGAL::Orthtree_traits_d` -*/ -template -struct Orthtree_traits_2 -{ -public: - - /// \name Types - /// @{ - - typedef Dimension_tag<2> Dimension; ///< Dimension type. - typedef Bbox_2 Bbox_d; ///< Bounding box type. - typedef typename GeomTraits::FT FT; ///< Number type. - typedef typename GeomTraits::Point_2 Point_d; ///< Point type. - typedef typename GeomTraits::Circle_2 Sphere_d; ///< Sphere type. - typedef typename GeomTraits::Cartesian_const_iterator_2 Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates. - typedef std::array Array; ///< Array type. - - /*! - * \brief Two directions along each axis in Cartesian space, relative to a node. - * - * Directions are mapped to numbers as 2-bit integers. - * - * The first bit indicates the axis (0 = x, 1 = y), - * the second bit indicates the direction along that axis (0 = -, 1 = +). - * - * The following diagram may be a useful reference: - * - * 3 * - * | - * | y+ - * | * - * 0 *------+------* 1 | - * | | - * | +-----* x+ - * | - * * 2 - * - * This lookup table may also be helpful: - * - * | Direction | bitset | number | Enum | - * | --------- | ------ | ------ | ----- | - * | `-x` | 00 | 0 | LEFT | - * | `+x` | 01 | 1 | RIGHT | - * | `-y` | 10 | 2 | DOWN | - * | `+y` | 11 | 3 | UP | - */ - enum Adjacency - { - LEFT, - RIGHT, - DOWN, - UP - }; - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Point_d` from an `Array` object. - */ - typedef unspecified_type Construct_point_d_from_array; -#else - struct Construct_point_d_from_array - { - Point_d operator() (const Array& array) const - { - return Point_d (array[0], array[1]); - } - }; -#endif - - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points). - */ - typedef unspecified_type Construct_bbox_d; -#else - struct Construct_bbox_d - { - Bbox_d operator() (const Array& min, - const Array& max) const - { - return Bbox_d (min[0], min[1], max[0], max[1]); - } - }; -#endif - - /// @} - - /// \name Operations - /// @{ - - /*! - Function used to construct an object of type `Construct_point_d_from_array`. - */ - Construct_point_d_from_array construct_point_d_from_array_object() const - { return Construct_point_d_from_array(); } - - /*! - Function used to construct an object of type `Construct_bbox_d`. - */ - Construct_bbox_d construct_bbox_d_object() const - { return Construct_bbox_d(); } - - /// @} -}; - -} - -#endif // CGAL_ORTHTREE_TRAITS_2_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_3.h b/Orthtree/include/CGAL/Orthtree_traits_3.h deleted file mode 100644 index c50028fb8f27..000000000000 --- a/Orthtree/include/CGAL/Orthtree_traits_3.h +++ /dev/null @@ -1,160 +0,0 @@ -// Copyright (c) 2020 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) : Simon Giraudot - -#ifndef CGAL_ORTHTREE_TRAITS_3_H -#define CGAL_ORTHTREE_TRAITS_3_H - -#include - -#include -#include - -namespace CGAL -{ - -/*! - \ingroup PkgOrthtreeTraits - - The class `Orthtree_traits_3` can be used as a template parameter of - the `Orthtree` class. - - \tparam GeomTraits model of `Kernel`. - - \cgalModels{OrthtreeTraits} - \sa `CGAL::Octree` - \sa `CGAL::Orthtree_traits_2` - \sa `CGAL::Orthtree_traits_d` -*/ -template -struct Orthtree_traits_3 -{ -public: - - /// \name Types - /// @{ - - typedef Dimension_tag<3> Dimension; ///< Dimension type. - typedef Bbox_3 Bbox_d; ///< Bounding box type. - typedef typename GeomTraits::FT FT; ///< Number type. - typedef typename GeomTraits::Point_3 Point_d; ///< Point type. - typedef typename GeomTraits::Sphere_3 Sphere_d; ///< Sphere type. - typedef typename GeomTraits::Cartesian_const_iterator_3 Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates. - typedef std::array Array; ///< Array type. - - /*! - * \brief Two directions along each axis in Cartesian space, relative to a node. - * - * Directions are mapped to numbers as 3-bit integers, - * though the numbers 6 and 7 are not used because there are only 6 different directions. - * - * The first two bits indicate the axis (00 = x, 01 = y, 10 = z), - * the third bit indicates the direction along that axis (0 = -, 1 = +). - * - * The following diagram may be a useful reference: - * - * 3 * - * | * 5 - * | / y+ - * |/ * z+ - * 0 *------+------* 1 | * - * /| |/ - * / | +-----* x+ - * 4 * | - * * 2 - * - * This lookup table may also be helpful: - * - * | Direction | bitset | number | Enum | - * | --------- | ------ | ------ | ----- | - * | `-x` | 000 | 0 | LEFT | - * | `+x` | 001 | 1 | RIGHT | - * | `-y` | 010 | 2 | DOWN | - * | `+y` | 011 | 3 | UP | - * | `-z` | 100 | 4 | BACK | - * | `+z` | 101 | 5 | FRONT | - */ - enum Adjacency - { - LEFT, - RIGHT, - DOWN, - UP, - BACK, - FRONT - }; - - /// \cond SKIP_IN_MANUAL - enum Child { - LEFT_BOTTOM_BACK, - RIGHT_BOTTOM_BACK, - LEFT_TOP_BACK, - RIGHT_TOP_BACK, - LEFT_BOTTOM_FRONT, - RIGHT_BOTTOM_FRONT, - LEFT_TOP_FRONT, - RIGHT_TOP_FRONT - }; - /// \endcond - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Point_d` from an `Array` object. - */ - typedef unspecified_type Construct_point_d_from_array; -#else - struct Construct_point_d_from_array - { - Point_d operator() (const Array& array) const - { - return Point_d (array[0], array[1], array[2]); - } - }; -#endif - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points). - */ - typedef unspecified_type Construct_bbox_d; -#else - struct Construct_bbox_d - { - Bbox_d operator() (const Array& min, - const Array& max) const - { - return Bbox_d (min[0], min[1], min[2], max[0], max[1], max[2]); - } - }; -#endif - - /// @} - - /// \name Operations - /// @{ - - /*! - Function used to construct an object of type `Construct_point_d_from_array`. - */ - Construct_point_d_from_array construct_point_d_from_array_object() const - { return Construct_point_d_from_array(); } - - /*! - Function used to construct an object of type `Construct_bbox_d`. - */ - Construct_bbox_d construct_bbox_d_object() const - { return Construct_bbox_d(); } - - /// @} -}; - -} - -#endif // CGAL_ORTHTREE_TRAITS_3_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_base.h b/Orthtree/include/CGAL/Orthtree_traits_base.h new file mode 100644 index 000000000000..9750b2684ac6 --- /dev/null +++ b/Orthtree/include/CGAL/Orthtree_traits_base.h @@ -0,0 +1,167 @@ +// Copyright (c) 2023 INRIA (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) : Jackson Campolattaro + +#ifndef ORTHTREE_ORTHTREE_TRAITS_BASE_H +#define ORTHTREE_ORTHTREE_TRAITS_BASE_H + +#include + +#include +#include +#include +#include + +namespace CGAL { + +/*! + \ingroup PkgOrthtreeTraits + + The class `Orthtree_traits_base` is a base class providing common choices for types and functors. + The base class is extended by `CGAL::Orthtree_traits_point` and by `CGAL::Orthtree_traits_face_graph`. + + \tparam GeomTraits a model of `Kernel`. + \tparam dim dimension of the ambient Euclidean space. + + \sa `CGAL::Orthtree_traits_point` + \sa `CGAL::Orthtree_traits_face_graph` +*/ + +template +struct Orthtree_traits_base { + /// \name Types + /// @{ + using Node_index = std::size_t; + using Kernel = GeomTraits; + static constexpr int dimension = dim; + using FT = typename GeomTraits::FT; + using Point_d = typename GeomTraits::Point_d; + using Bbox_d = typename GeomTraits::Iso_box_d; + using Sphere_d = typename GeomTraits::Sphere_d; + using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_d; + /*! + * Adjacency type. + * + * \note This type is used to identify adjacency directions with + * easily understandable keywords (left, right, up, down, ...) and is thus + * mainly useful in 2D and 3D. In + * higher dimensions, such keywords do not exist and this type is + * simply an integer. Conversions from this integer to bitsets still + * work but do not provide any user-friendly API for adjacency selection. + * + * Two directions along each axis in %Cartesian space, relative to a node. + * + * Directions are mapped to numbers as 3-bit integers in the 3D case or as 2-bit integers in the 2D case. + * In the 3d case the numbers 6 and 7 are not used because there are only 6 different directions. + * + * The first two bits indicate the axis (00 = x, 01 = y, 10 = z), + * the third bit indicates the direction along that axis (0 = -, 1 = +). + * + * The following diagram and table showing the 3D case may be a useful reference (2D case is identical with one dimension less): + * + * 3 * + * | * 4 + * | / y+ + * |/ * + * 0 *------+------* 1 | + * /| | + * / | +-----* x+ + * 5 * | / + * * 2 / + * * z+ + * + * This lookup table may also be helpful: + * + * | Direction | bitset | number | Enum | + * | --------- | ------ | ------ | ----- | + * | `-x` | 000 | 0 | LEFT | + * | `+x` | 001 | 1 | RIGHT | + * | `-y` | 010 | 2 | DOWN | + * | `+y` | 011 | 3 | UP | + * | `-z` | 100 | 4 | BACK | + * | `+z` | 101 | 5 | FRONT | + */ + using Adjacency = int; + /// @} + + auto construct_point_d_object() const { + return [](auto... Args) -> Point_d { + std::initializer_list args_list{Args...}; + return Point_d{static_cast(args_list.size()), args_list.begin(), args_list.end()}; + }; + } +}; + +template +struct Orthtree_traits_base { + using Node_index = std::size_t; + using Kernel = GeomTraits; + static constexpr int dimension = 2; + using FT = typename GeomTraits::FT; + using Point_d = typename GeomTraits::Point_2; + using Bbox_d = typename GeomTraits::Iso_rectangle_2; + using Sphere_d = typename GeomTraits::Circle_2; + using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_2; + + enum Adjacency { + LEFT, + RIGHT, + DOWN, + UP + }; + + auto construct_point_d_object() const { + return [](const FT& x, const FT& y) -> Point_d { + return {x, y}; + }; + } +}; + +template +struct Orthtree_traits_base { + using Node_index = std::size_t; + using Kernel = GeomTraits; + static constexpr int dimension = 3; + using FT = typename GeomTraits::FT; + using Point_d = typename GeomTraits::Point_3; + using Bbox_d = typename GeomTraits::Iso_cuboid_3; + using Sphere_d = typename GeomTraits::Sphere_3; + using Cartesian_const_iterator_d = typename GeomTraits::Cartesian_const_iterator_3; + + enum Adjacency { + LEFT, + RIGHT, + DOWN, + UP, + BACK, + FRONT + }; + + enum Child { + LEFT_BOTTOM_BACK, + RIGHT_BOTTOM_BACK, + LEFT_TOP_BACK, + RIGHT_TOP_BACK, + LEFT_BOTTOM_FRONT, + RIGHT_BOTTOM_FRONT, + LEFT_TOP_FRONT, + RIGHT_TOP_FRONT + }; + + auto construct_point_d_object() const { + return [](const FT& x, const FT& y, const FT& z) -> Point_d { + return {x, y, z}; + }; + } +}; + +} + +#endif //ORTHTREE_ORTHTREE_TRAITS_BASE_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_d.h b/Orthtree/include/CGAL/Orthtree_traits_d.h deleted file mode 100644 index 4f48caba07ff..000000000000 --- a/Orthtree/include/CGAL/Orthtree_traits_d.h +++ /dev/null @@ -1,137 +0,0 @@ -// Copyright (c) 2020 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) : Simon Giraudot - -#ifndef CGAL_ORTHTREE_TRAITS_D_H -#define CGAL_ORTHTREE_TRAITS_D_H - -#include - -#include - -namespace CGAL -{ - -/*! - \ingroup PkgOrthtreeTraits - - The class `Orthtree_traits_d` can be used as a template parameter of - the `Orthtree` class. - - \tparam GeomTraits model of `Kernel`. - \tparam DimensionTag specialization of `CGAL::Dimension_tag`. - - \cgalModels{OrthtreeTraits} - \sa `CGAL::Orthtree` - \sa `CGAL::Orthtree_traits_2` - \sa `CGAL::Orthtree_traits_3` -*/ - -template -struct Orthtree_traits_d -{ -public: - - /// \name Types - /// @{ - - typedef DimensionTag Dimension; ///< Dimension type. - typedef typename GeomTraits::FT FT; ///< Number type. - typedef typename GeomTraits::Point_d Point_d; ///< Point type. - typedef typename GeomTraits::Sphere_d Sphere_d; ///< Sphere type. - typedef typename GeomTraits::Cartesian_const_iterator_d Cartesian_const_iterator_d; ///< An iterator over the %Cartesian coordinates. - typedef std::array Array; ///< Array type. - -#ifdef DOXYGEN_RUNNING - typedef unspecified_type Bbox_d; ///< Bounding box type. -#else - class Bbox_d - { - Point_d m_min, m_max; - public: - - Bbox_d (const Point_d& pmin, const Point_d& pmax) - : m_min (pmin), m_max (pmax) - { } - - const Point_d& min BOOST_PREVENT_MACRO_SUBSTITUTION () { return m_min; } - const Point_d& max BOOST_PREVENT_MACRO_SUBSTITUTION () { return m_max; } - }; -#endif - - /*! - Adjacency type. - - \note This type is used to identify adjacency directions with - easily understandable keywords (left, right, up, etc.) and is thus - mainly useful for `Orthtree_traits_2` and `Orthtree_traits_3`. In - higher dimensions, such keywords do not exist and this type is - simply an integer. Conversions from this integer to bitsets still - works but do not provide any easier API for adjacency selection. - */ - typedef int Adjacency; - - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Point_d` from an `Array` object. - */ - typedef unspecified_type Construct_point_d_from_array; -#else - struct Construct_point_d_from_array - { - Point_d operator() (const Array& array) const - { - return Point_d (array.begin(), array.end()); - } - }; -#endif - - -#ifdef DOXYGEN_RUNNING - /*! - Functor with an operator to construct a `Bbox_d` from two `Array` objects (coordinates of minimum and maximum points). - */ - typedef unspecified_type Construct_bbox_d; -#else - struct Construct_bbox_d - { - Bbox_d operator() (const Array& min, - const Array& max) const - { - return Bbox_d (Point_d (min.begin(), min.end()), - Point_d (max.begin(), max.end())); - } - }; -#endif - - /// @} - - /// \name Operations - /// @{ - - /*! - Function used to construct an object of type `Construct_point_d_from_array`. - */ - Construct_point_d_from_array construct_point_d_from_array_object() const - { return Construct_point_d_from_array(); } - - /*! - Function used to construct an object of type `Construct_bbox_d`. - */ - Construct_bbox_d construct_bbox_d_object() const - { return Construct_bbox_d(); } - - /// @} -}; - -} - -#endif // CGAL_ORTHTREE_TRAITS_D_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_face_graph.h b/Orthtree/include/CGAL/Orthtree_traits_face_graph.h new file mode 100644 index 000000000000..ffa353cd6133 --- /dev/null +++ b/Orthtree/include/CGAL/Orthtree_traits_face_graph.h @@ -0,0 +1,168 @@ +// Copyright (c) 2023 INRIA (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) : Sebastien LORIOT + + +#ifndef CGAL_ORTHREE_TRAITS_FACE_GRAPH_H +#define CGAL_ORTHREE_TRAITS_FACE_GRAPH_H + +#include + +#include + +#include +#include +#include + +namespace CGAL { + +/*! +\ingroup PkgOrthtreeTraits + +Traits class for the `Orthtree` class to be used to construct a 3D octree around +a triangulated surface mesh. Each node of the octree will store all the faces of the +mesh intersected by its bounding box. The subdivision of the octree is controlled +by the nested class `Orthtree_traits_face_graph::Split_predicate_node_min_extent` +to which the minimal extent of a node should be provided. + +\tparam TriangleMesh a model of `FaceListGraph` with all faces being triangles +\tparam VertexPointMap a property map associating points to the vertices of `TriangleMesh` + +\cgalModels{OrthtreeTraitsWithData} +*/ +template +struct Orthtree_traits_face_graph : public Orthtree_traits_base< + typename Kernel_traits::value_type>::type, 3 > { + + Orthtree_traits_face_graph(const TriangleMesh& pm, VertexPointMap vpm) + : m_pm(pm), m_vpm(vpm) {} + + /// \name Types + /// @{ + + using Base = Orthtree_traits_base< + typename Kernel_traits::value_type>::type, 3 >; + using Self = Orthtree_traits_face_graph; + using Tree = Orthtree; + + using Point_d = typename Self::Point_d; + using Bbox_d = typename Self::Bbox_d; + using FT = typename Self::FT; + using Cartesian_const_iterator_d = typename Self::Cartesian_const_iterator_d; + + using Node_index = typename Base::Node_index; + using Node_data = std::vector::face_descriptor>; + + using Geom_traits = typename Kernel_traits::type; + + using Construct_root_node_bbox = std::function; + using Construct_root_node_contents = std::function; + using Distribute_node_contents = std::function; + + /// @} + + /// \name Operations + /// @{ + + auto construct_root_node_bbox_object() const { + return [&]() -> Bbox_d { + + std::array min = {0.0, 0}, max = {0.0, 0}; + if (faces(m_pm).begin() != faces(m_pm).end()) { + bool first = true; + for (auto v: vertices(m_pm)) { + const Point_d& p_v = get(m_vpm, v); + for (int i = 0; i < 3; ++i) { + if (first || p_v[i] < min[i]) min[i] = p_v[i]; + if (first || p_v[i] > max[i]) max[i] = p_v[i]; + } + first=false; + } + } + + return {std::apply(Self::construct_point_d_object(), min), + std::apply(Self::construct_point_d_object(), max)}; + }; + } + + auto construct_root_node_contents_object() const { + return [&]() -> Node_data { + return {faces(m_pm).begin(), faces(m_pm).end()}; + }; + } + + auto distribute_node_contents_object() const { + return [&](Node_index n, Tree& tree, const Point_d& /* center */) -> void { + Node_data& ndata = tree.data(n); + for (int i = 0; i < 8; ++i) { + Node_index child = tree.child(n, i); + Node_data& child_data = tree.data(child); + Bbox_d bbox = tree.bbox(child); + for (auto f : ndata) { + typename boost::graph_traits::halfedge_descriptor + h = halfedge(f, m_pm); + typename Geom_traits::Triangle_3 t(get(m_vpm, source(h, m_pm)), + get(m_vpm, target(h, m_pm)), + get(m_vpm, target(next(h, m_pm), m_pm))); + if (do_intersect(t, bbox)) + child_data.push_back(f); + } + } + }; + } + + /// @} + + /// Recommended split predicate to pass to `Orthtree::refine()` function so + /// that the octree is refined until a node is either empty or has an extent + /// that would be smaller after split than the corresponding value provided to the constructor. + class Split_predicate_node_min_extent { + + std::array m_min_extent; + + public: + + /// constructor with `me` being the minimal value a node extent could be + /// (same value for all dimension). + Split_predicate_node_min_extent(const FT& me) + : m_min_extent({me, me, me}) {} + + /// constructor with `me` being the minimal value a node extent could be + /// (one value per dimension). + Split_predicate_node_min_extent(const std::array& me) + : m_min_extent(me) {} + + /*! + \brief returns `true` if `ni` should be split, `false` otherwise. + */ + template + bool operator()(NodeIndex ni, const Tree& tree) const { + if (tree.data(ni).empty()) return false; + + Bbox_d bb = tree.bbox(ni); + + for (int i = 0; i < 3; ++i) + if (((bb.max)()[i] - (bb.min)()[i]) < 2 * m_min_extent[i]) + return false; + return true; + } + }; + + +private: + + const TriangleMesh& m_pm; + VertexPointMap m_vpm; +}; + +} // end of CGAL namespace + + +#endif // CGAL_ORTHREE_TRAITS_FACE_GRAPH_H diff --git a/Orthtree/include/CGAL/Orthtree_traits_point.h b/Orthtree/include/CGAL/Orthtree_traits_point.h new file mode 100644 index 000000000000..84e95bde81c0 --- /dev/null +++ b/Orthtree/include/CGAL/Orthtree_traits_point.h @@ -0,0 +1,201 @@ +// Copyright (c) 2023 INRIA (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) : Jackson Campolattaro + +#ifndef ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H +#define ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H + +#include + +#include +#include +#include +#include + +#include + +namespace CGAL { + +template +void reassign_points( + Tree& tree, PointMap& point_map, + typename Tree::Node_index n, const typename Tree::Point& center, typename Tree::Node_data points, + std::bitset coord = {}, std::size_t dimension = 0 +) { + + // Root case: reached the last dimension + if (dimension == Tree::dimension) { + tree.data(tree.child(n, coord.to_ulong())) = points; + return; + } + + // Split the point collection around the center point on this dimension + auto split_point = std::partition( + points.begin(), points.end(), + [&](const auto& p) -> bool { + return (get(point_map, p)[int(dimension)] < center[int(dimension)]); + } + ); + + // Further subdivide the first side of the split + std::bitset coord_left = coord; + coord_left[dimension] = false; + reassign_points(tree, point_map, n, center, {points.begin(), split_point}, coord_left, dimension + 1); + + // Further subdivide the second side of the split + std::bitset coord_right = coord; + coord_right[dimension] = true; + reassign_points(tree, point_map, n, center, {split_point, points.end()}, coord_right, dimension + 1); +} + +/*! + \ingroup PkgOrthtreeTraits + + Traits class for defining an orthtree of points using the class `CGAL::Orthtree`. + + \tparam GeomTraits model of `Kernel`. + \tparam PointRange must be a model of `Range` whose value type is the key type of `PointMap` and whose iterator type is model of `RandomAccessIterator` + \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is a point type from `GeomTraits` matching the current dimension + \tparam dimension the dimension of the ambient Euclidean space. + + \warning The input point set is not copied. It is used directly + and is rearranged by the `Orthtree`. Altering the point range + after creating the orthtree will leave it in an invalid state. + + \cgalModels{CollectionPartitioningOrthtreeTraits} + \sa `CGAL::Octree` + \sa `CGAL::Quadtree` + \sa `CGAL::Orthtree_traits_base` +*/ +template < + typename GeomTraits, + typename PointRange, + typename PointMap = Identity_property_map::value_type>, + bool hypercubic_nodes = false, + int dimension = Ambient_dimension< + typename std::iterator_traits::value_type, + GeomTraits + >::value +> +struct Orthtree_traits_point : public Orthtree_traits_base { +public: + /// \name Types + /// @{ + using Node_data = boost::iterator_range; + /// @} + + using Base = Orthtree_traits_base; + using Self = Orthtree_traits_point; + using Tree = Orthtree; + + using Node_index = typename Base::Node_index; + using Node_data_element = typename std::iterator_traits::value_type; + + Orthtree_traits_point( + PointRange& points, + PointMap point_map = PointMap() + ) : m_points(points), m_point_map(point_map) {} + + auto construct_root_node_bbox_object() const { + return [&]() -> typename Self::Bbox_d { + + std::array bbox_min, bbox_max; + Orthtrees::internal::Cartesian_ranges cartesian_range; + + // init bbox with first values found + { + const typename Self::Point_d& point = get(m_point_map, *(m_points.begin())); + std::size_t i = 0; + for (const typename Self::FT& x: cartesian_range(point)) { + bbox_min[i] = x; + bbox_max[i] = x; + ++i; + } + } + // Expand bbox to contain all points + for (const auto& p: m_points) { + const typename Self::Point_d& point = get(m_point_map, p); + std::size_t i = 0; + for (const typename Self::FT& x: cartesian_range(point)) { + bbox_min[i] = (std::min)(x, bbox_min[i]); + bbox_max[i] = (std::max)(x, bbox_max[i]); + ++i; + } + } + +#if !defined(_MSC_VER) || _MSC_VER > 1920 + if constexpr (hypercubic_nodes) { +#else + if (hypercubic_nodes) { +#endif + std::array center; + typename Self::FT max_side = 0; + for (int i = 0; i < Self::dimension; i++) { + typename Self::FT side = bbox_max[i] - bbox_min[i]; + max_side = (std::max)(max_side, side); + center[i] = (bbox_min[i] + bbox_max[i]) * 0.5f; + } + max_side *= 0.5f; + for (int i = 0; i < Self::dimension; i++) { + bbox_min[i] = center[i] - max_side; + bbox_max[i] = center[i] + max_side; + } + } + + return {std::apply(Self::construct_point_d_object(), bbox_min), + std::apply(Self::construct_point_d_object(), bbox_max)}; + }; + } + + auto construct_root_node_contents_object() const { + return [&]() -> typename Self::Node_data { + return {m_points.begin(), m_points.end()}; + }; + } + + auto distribute_node_contents_object() const { + return [&](Node_index n, Tree& tree, const typename Self::Point_d& center) { + CGAL_precondition(!tree.is_leaf(n)); + reassign_points(tree, m_point_map, n, center, tree.data(n)); + }; + } + + auto construct_sphere_d_object() const { + return [](const typename Self::Point_d& center, const typename Self::FT& squared_radius) -> typename Self::Sphere_d { + return typename Self::Sphere_d(center, squared_radius); + }; + } + + auto construct_center_d_object() const { + return [](const typename Self::Sphere_d& sphere) -> typename Self::Point_d { + return sphere.center(); + }; + } + + auto compute_squared_radius_d_object() const { + return [](const typename Self::Sphere_d& sphere) -> typename Self::FT { + return sphere.squared_radius(); + }; + } + + auto squared_distance_of_element_object() const { + return [&](const Node_data_element& index, const typename Self::Point_d& point) -> typename Self::FT { + return CGAL::squared_distance(get(m_point_map, index), point); + }; + } + + PointRange& m_points; + PointMap m_point_map; +}; + +} + + +#endif //ORTHTREE_TESTS_ORTHTREE_TRAITS_POINT_H diff --git a/Orthtree/include/CGAL/Quadtree.h b/Orthtree/include/CGAL/Quadtree.h index 0434df2d13ab..5b6086bcf73e 100644 --- a/Orthtree/include/CGAL/Quadtree.h +++ b/Orthtree/include/CGAL/Quadtree.h @@ -15,34 +15,28 @@ #include #include -#include +#include namespace CGAL { /*! - \ingroup PkgOrthtreeClasses + \ingroup PkgOrthtreeRef - \brief Alias that specializes the `Orthtree` class to a 2D quadtree. - - These two types are exactly equivalent: - - `Quadtree` - - `Orthtree, PointRange, PointMap>`. - - \warning This is a not a real class but an alias, please refer to - the documentation of `Orthtree`. + \brief Alias that specializes the `Orthtree` class to a 2D quadtree storing 2D points. \tparam GeomTraits must be a model of `Kernel` - \tparam PointRange_ must be a model of range whose value type is the key type of `PointMap` + \tparam PointRange must be a model of `Range` whose value type is the key type of `PointMap` \tparam PointMap must be a model of `ReadablePropertyMap` whose value type is `GeomTraits::Point_2` + \tparam square_nodes Boolean to enforce square nodes */ template ::value_type> > -#ifdef DOXYGEN_RUNNING -class Quadtree; -#else -using Quadtree = Orthtree, PointRange, PointMap>; -#endif + ::value_type>, + bool squared_nodes = false +> + +using Quadtree = Orthtree>; + } // namespace CGAL diff --git a/Orthtree/package_info/Orthtree/dependencies b/Orthtree/package_info/Orthtree/dependencies index bcb222ae1968..d8fbb010a7a2 100644 --- a/Orthtree/package_info/Orthtree/dependencies +++ b/Orthtree/package_info/Orthtree/dependencies @@ -1,6 +1,10 @@ Algebraic_foundations +Cartesian_kernel +Circulator Distance_2 Distance_3 +Filtered_kernel +Hash_map Installation Intersections_2 Intersections_3 @@ -9,7 +13,11 @@ Kernel_23 Modular_arithmetic Number_types Orthtree +Point_set_2 Profiling_tools Property_map STL_Extension +Spatial_sorting Stream_support +TDS_2 +Triangulation_2 diff --git a/Orthtree/test/Orthtree/CMakeLists.txt b/Orthtree/test/Orthtree/CMakeLists.txt index 6776af021578..c483d861c739 100644 --- a/Orthtree/test/Orthtree/CMakeLists.txt +++ b/Orthtree/test/Orthtree/CMakeLists.txt @@ -16,6 +16,7 @@ create_single_source_cgal_program("test_octree_traverse.cpp") create_single_source_cgal_program("test_octree_intersecting.cpp") create_single_source_cgal_program("test_octree_copy_move_constructors.cpp") create_single_source_cgal_program("test_octree_kernels.cpp") +create_single_source_cgal_program("test_octree_custom_properties.cpp") create_single_source_cgal_program("test_node_index.cpp") create_single_source_cgal_program("test_node_adjacent.cpp") diff --git a/Orthtree/test/Orthtree/test_node.cpp b/Orthtree/test/Orthtree/test_node.cpp deleted file mode 100644 index 13e3a3f07285..000000000000 --- a/Orthtree/test/Orthtree/test_node.cpp +++ /dev/null @@ -1,55 +0,0 @@ - -#include -#include -#include -#include - -typedef CGAL::Orthtree::Node::iterator> Node; - -int main(void) { - - // Build a new node - Node n = Node(); - - // Check that its values are correct - assert(n.is_root()); - assert(n.is_leaf()); - assert(!n.parent()); - assert(n.depth() == 0); - assert(n.location()[0] == 0 && n.location()[1] == 0 && n.location()[2] == 0); - - // Split the node - n.split(); - - // Check that it's children's values are also correct - for (std::size_t i = 0; i < 8; ++i) { - - assert(!n[i].is_root()); - assert(n[i].is_leaf()); - assert(*n[i].parent() == n); - assert(n[i].depth() == 1); - } - - // Check that the parent has updated - assert(n.is_root()); - assert(!n.is_leaf()); - - // Split one of the children - n[1].split(); - - // Check each of that child's children - for (std::size_t i = 0; i < 8; ++i) { - - assert(!n[1][i].is_root()); - assert(n[1][i].is_leaf()); - assert(*n[1][i].parent() == n[1]); - assert(*n[1][i].parent()->parent() == n); - assert(n[1][i].depth() == 2); - } - - // Check that the child's values have updated - assert(!n[1].is_root()); - assert(!n[1].is_leaf()); - - return 0; -} diff --git a/Orthtree/test/Orthtree/test_node_adjacent.cpp b/Orthtree/test/Orthtree/test_node_adjacent.cpp index 084d72dc5176..db9581b7225b 100644 --- a/Orthtree/test/Orthtree/test_node_adjacent.cpp +++ b/Orthtree/test/Orthtree/test_node_adjacent.cpp @@ -3,16 +3,14 @@ #include #include -#include #include +#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; -typedef Octree::Node Node; -typedef Octree::Traits Traits; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; +using Traits = Octree::Traits; int main(void) { @@ -42,33 +40,43 @@ int main(void) { std::cout << octree << std::endl; // Root node should have no siblings - assert(octree.root().adjacent_node(0).is_null()); - assert(octree.root().adjacent_node(1).is_null()); - assert(octree.root().adjacent_node(2).is_null()); - assert(octree.root().adjacent_node(3).is_null()); - assert(octree.root().adjacent_node(4).is_null()); - assert(octree.root().adjacent_node(5).is_null()); + assert(!octree.adjacent_node(octree.root(), 0)); + assert(!octree.adjacent_node(octree.root(), 1)); + assert(!octree.adjacent_node(octree.root(), 2)); + assert(!octree.adjacent_node(octree.root(), 3)); + assert(!octree.adjacent_node(octree.root(), 4)); + assert(!octree.adjacent_node(octree.root(), 5)); // Left Top Front node should have siblings to the Right, Down, and Back - auto left_top_back = octree.root()[Traits::LEFT_TOP_BACK]; + auto left_top_back = octree.node(Traits::LEFT_TOP_BACK); + + assert(octree.node(Traits::RIGHT_TOP_BACK) == + *octree.adjacent_node(left_top_back, Traits::RIGHT)); + assert(octree.node(Traits::LEFT_BOTTOM_BACK) == + *octree.adjacent_node(left_top_back, Traits::DOWN)); + assert(octree.node(Traits::LEFT_TOP_FRONT) == + *octree.adjacent_node(left_top_back, Traits::FRONT)); + assert(!octree.adjacent_node(left_top_back, Traits::LEFT)); + assert(!octree.adjacent_node(left_top_back, Traits::UP)); + assert(!octree.adjacent_node(left_top_back, Traits::BACK)); - assert(octree.root()[Traits::RIGHT_TOP_BACK] == left_top_back.adjacent_node(Traits::RIGHT)); - assert(octree.root()[Traits::LEFT_BOTTOM_BACK] == left_top_back.adjacent_node(Traits::DOWN)); - assert(octree.root()[Traits::LEFT_TOP_FRONT] == left_top_back.adjacent_node(Traits::FRONT)); - assert(left_top_back.adjacent_node(Traits::LEFT).is_null()); - assert(left_top_back.adjacent_node(Traits::UP).is_null()); - assert(left_top_back.adjacent_node(Traits::BACK).is_null()); + auto right_top_back_of_left_bottom_back = octree.node(Traits::LEFT_BOTTOM_BACK, Traits::RIGHT_TOP_BACK); - std::cout << octree.root()[Traits::LEFT_BOTTOM_BACK] << std::endl; + assert( + octree.node(Traits::LEFT_BOTTOM_BACK, Traits::LEFT_TOP_BACK) == + octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::LEFT) + ); + assert( + octree.node(Traits::RIGHT_BOTTOM_BACK) == + octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::RIGHT) + ); + assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::RIGHT).has_value()); + assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::UP).has_value()); + assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::DOWN).has_value()); + assert(octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::FRONT).has_value()); - auto right_top_back_of_left_bottom_back = octree.root()[Traits::LEFT_BOTTOM_BACK][Traits::RIGHT_TOP_BACK]; - assert(octree.root()[Traits::LEFT_BOTTOM_BACK][Traits::LEFT_TOP_BACK] == right_top_back_of_left_bottom_back.adjacent_node(Traits::LEFT)); - assert(octree.root()[Traits::RIGHT_BOTTOM_BACK] == right_top_back_of_left_bottom_back.adjacent_node(Traits::RIGHT)); - assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::RIGHT).is_null()); - assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::UP).is_null()); - assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::DOWN).is_null()); - assert(right_top_back_of_left_bottom_back.adjacent_node(Traits::BACK).is_null()); - assert(!right_top_back_of_left_bottom_back.adjacent_node(Traits::FRONT).is_null()); + // A node at the back of the tree should have no neighbor to its back + assert(!octree.adjacent_node(right_top_back_of_left_bottom_back, Traits::BACK)); return 0; } diff --git a/Orthtree/test/Orthtree/test_node_index.cpp b/Orthtree/test/Orthtree/test_node_index.cpp index 9ee08bd1083c..d8216f6e996a 100644 --- a/Orthtree/test/Orthtree/test_node_index.cpp +++ b/Orthtree/test/Orthtree/test_node_index.cpp @@ -3,14 +3,13 @@ #include #include -#include #include +#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; int main(void) { @@ -37,12 +36,11 @@ int main(void) { Octree octree(points, points.point_map()); octree.refine(10, 1); - std::cout << "root: " << octree.root().local_coordinates() << std::endl; - std::cout << "first child: " << octree.root()[0].local_coordinates() << std::endl; - std::cout << "fifth child: " << octree.root()[4].local_coordinates() << std::endl; - std::cout << "fifth child of first child: " << octree.root()[0][4].local_coordinates() << std::endl; - - // TODO + std::cout << "root: " << octree.local_coordinates(octree.root()) << std::endl; + std::cout << "first child: " << octree.local_coordinates(octree.child(octree.root(), 0)) << std::endl; + std::cout << "fifth child: " << octree.local_coordinates(octree.child(octree.root(), 4)) << std::endl; + std::cout << "fifth child of first child: " + << octree.local_coordinates(octree.child(octree.child(octree.root(), 0), 4)) << std::endl; return 0; } diff --git a/Orthtree/test/Orthtree/test_octree_bbox.cpp b/Orthtree/test/Orthtree/test_octree_bbox.cpp index 375894cc54cc..f36ba6aeb942 100644 --- a/Orthtree/test/Orthtree/test_octree_bbox.cpp +++ b/Orthtree/test/Orthtree/test_octree_bbox.cpp @@ -1,17 +1,14 @@ - -#include #include -#include #include - +#include +#include #include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; void test_1_node() { @@ -23,8 +20,10 @@ void test_1_node() { Octree octree(points, points.point_map()); octree.refine(10, 1); + Octree::Bbox expected_bbox{-1, -1, -1, -1, -1, -1}; + // Compare the top (only) node - assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1, -1, -1, -1, -1, -1)); + assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, -1, -1, -1)); } void test_9_nodes() { @@ -35,21 +34,21 @@ void test_9_nodes() { points.insert({1, 1, 1}); // Create the octree - Octree octree(points, points.point_map(), 1.1); + Octree octree(points, points.point_map()); octree.refine(10, 1); // Compare the top node - assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1.1, -1.1, -1.1, 1.1, 1.1, 1.1)); + assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, 1, 1, 1)); // Compare the child nodes - assert(octree.bbox(octree.root()[0]) == CGAL::Bbox_3(-1.1, -1.1, -1.1, 0, 0, 0)); - assert(octree.bbox(octree.root()[1]) == CGAL::Bbox_3(0, -1.1, -1.1, 1.1, 0, 0)); - assert(octree.bbox(octree.root()[2]) == CGAL::Bbox_3(-1.1, 0, -1.1, 0, 1.1, 0)); - assert(octree.bbox(octree.root()[3]) == CGAL::Bbox_3(0, 0, -1.1, 1.1, 1.1, 0)); - assert(octree.bbox(octree.root()[4]) == CGAL::Bbox_3(-1.1, -1.1, 0, 0, 0, 1.1)); - assert(octree.bbox(octree.root()[5]) == CGAL::Bbox_3(0, -1.1, 0, 1.1, 0, 1.1)); - assert(octree.bbox(octree.root()[6]) == CGAL::Bbox_3(-1.1, 0, 0, 0, 1.1, 1.1)); - assert(octree.bbox(octree.root()[7]) == CGAL::Bbox_3(0, 0, 0, 1.1, 1.1, 1.1)); + assert(octree.bbox(octree.node(0)) == Octree::Bbox(-1, -1, -1, 0, 0, 0)); + assert(octree.bbox(octree.node(1)) == Octree::Bbox(0, -1, -1, 1, 0, 0)); + assert(octree.bbox(octree.node(2)) == Octree::Bbox(-1, 0, -1, 0, 1, 0)); + assert(octree.bbox(octree.node(3)) == Octree::Bbox(0, 0, -1, 1, 1, 0)); + assert(octree.bbox(octree.node(4)) == Octree::Bbox(-1, -1, 0, 0, 0, 1)); + assert(octree.bbox(octree.node(5)) == Octree::Bbox(0, -1, 0, 1, 0, 1)); + assert(octree.bbox(octree.node(6)) == Octree::Bbox(-1, 0, 0, 0, 1, 1)); + assert(octree.bbox(octree.node(7)) == Octree::Bbox(0, 0, 0, 1, 1, 1)); } void test_25_nodes() { @@ -62,41 +61,70 @@ void test_25_nodes() { points.insert({1, 0.5, 1}); // Create the octree - Octree octree(points, points.point_map(), 1.5); + Octree octree(points, points.point_map()); octree.refine(10, 1); // Compare the top node - assert(octree.bbox(octree.root()) == CGAL::Bbox_3(-1.5, -1.5, -1.5, 1.5, 1.5, 1.5)); + assert(octree.bbox(octree.root()) == Octree::Bbox(-1, -1, -1, 1, 1, 1)); // Compare the child nodes - assert(octree.bbox(octree.root()[0]) == CGAL::Bbox_3(-1.5, -1.5, -1.5, 0, 0, 0)); - assert(octree.bbox(octree.root()[1]) == CGAL::Bbox_3(0, -1.5, -1.5, 1.5, 0, 0)); - assert(octree.bbox(octree.root()[2]) == CGAL::Bbox_3(-1.5, 0, -1.5, 0, 1.5, 0)); - assert(octree.bbox(octree.root()[3]) == CGAL::Bbox_3(0, 0, -1.5, 1.5, 1.5, 0)); - assert(octree.bbox(octree.root()[4]) == CGAL::Bbox_3(-1.5, -1.5, 0, 0, 0, 1.5)); - assert(octree.bbox(octree.root()[5]) == CGAL::Bbox_3(0, -1.5, 0, 1.5, 0, 1.5)); - assert(octree.bbox(octree.root()[6]) == CGAL::Bbox_3(-1.5, 0, 0, 0, 1.5, 1.5)); - assert(octree.bbox(octree.root()[7]) == CGAL::Bbox_3(0, 0, 0, 1.5, 1.5, 1.5)); + assert(octree.bbox(octree.node(0)) == Octree::Bbox(-1, -1, -1, 0, 0, 0)); + assert(octree.bbox(octree.node(1)) == Octree::Bbox(0, -1, -1, 1, 0, 0)); + assert(octree.bbox(octree.node(2)) == Octree::Bbox(-1, 0, -1, 0, 1, 0)); + assert(octree.bbox(octree.node(3)) == Octree::Bbox(0, 0, -1, 1, 1, 0)); + assert(octree.bbox(octree.node(4)) == Octree::Bbox(-1, -1, 0, 0, 0, 1)); + assert(octree.bbox(octree.node(5)) == Octree::Bbox(0, -1, 0, 1, 0, 1)); + assert(octree.bbox(octree.node(6)) == Octree::Bbox(-1, 0, 0, 0, 1, 1)); + assert(octree.bbox(octree.node(7)) == Octree::Bbox(0, 0, 0, 1, 1, 1)); // Compare children of the first child - assert(octree.bbox(octree.root()[0][0]) == CGAL::Bbox_3(-1.5, -1.5, -1.5, -0.75, -0.75, -0.75)); - assert(octree.bbox(octree.root()[0][1]) == CGAL::Bbox_3(-0.75, -1.5, -1.5, 0, -0.75, -0.75)); - assert(octree.bbox(octree.root()[0][2]) == CGAL::Bbox_3(-1.5, -0.75, -1.5, -0.75, 0, -0.75)); - assert(octree.bbox(octree.root()[0][3]) == CGAL::Bbox_3(-0.75, -0.75, -1.5, 0, 0, -0.75)); - assert(octree.bbox(octree.root()[0][4]) == CGAL::Bbox_3(-1.5, -1.5, -0.75, -0.75, -0.75, 0)); - assert(octree.bbox(octree.root()[0][5]) == CGAL::Bbox_3(-0.75, -1.5, -0.75, 0, -0.75, 0)); - assert(octree.bbox(octree.root()[0][6]) == CGAL::Bbox_3(-1.5, -0.75, -0.75, -0.75, 0, 0)); - assert(octree.bbox(octree.root()[0][7]) == CGAL::Bbox_3(-0.75, -0.75, -0.75, 0, 0, 0)); + assert(octree.bbox(octree.node(0, 0)) == + Octree::Bbox(-1, -1, -1, -0.5, -0.5, -0.5)); + assert(octree.bbox(octree.node(0, 1)) == + Octree::Bbox(-0.5, -1, -1, 0, -0.5, -0.5)); + assert(octree.bbox(octree.node(0, 2)) == + Octree::Bbox(-1, -0.5, -1, -0.5, 0, -0.5)); + assert(octree.bbox(octree.node(0, 3)) == + Octree::Bbox(-0.5, -0.5, -1, 0, 0, -0.5)); + assert(octree.bbox(octree.node(0, 4)) == + Octree::Bbox(-1, -1, -0.5, -0.5, -0.5, 0)); + assert(octree.bbox(octree.node(0, 5)) == + Octree::Bbox(-0.5, -1, -0.5, 0, -0.5, 0)); + assert(octree.bbox(octree.node(0, 6)) == + Octree::Bbox(-1, -0.5, -0.5, -0.5, 0, 0)); + assert(octree.bbox(octree.node(0, 7)) == + Octree::Bbox(-0.5, -0.5, -0.5, 0, 0, 0)); // Compare children of the last child - assert(octree.bbox(octree.root()[7][0]) == CGAL::Bbox_3(0, 0, 0, 0.75, 0.75, 0.75)); - assert(octree.bbox(octree.root()[7][1]) == CGAL::Bbox_3(0.75, 0, 0, 1.5, 0.75, 0.75)); - assert(octree.bbox(octree.root()[7][2]) == CGAL::Bbox_3(0, 0.75, 0, 0.75, 1.5, 0.75)); - assert(octree.bbox(octree.root()[7][3]) == CGAL::Bbox_3(0.75, 0.75, 0, 1.5, 1.5, 0.75)); - assert(octree.bbox(octree.root()[7][4]) == CGAL::Bbox_3(0, 0, 0.75, 0.75, 0.75, 1.5)); - assert(octree.bbox(octree.root()[7][5]) == CGAL::Bbox_3(0.75, 0, 0.75, 1.5, 0.75, 1.5)); - assert(octree.bbox(octree.root()[7][6]) == CGAL::Bbox_3(0, 0.75, 0.75, 0.75, 1.5, 1.5)); - assert(octree.bbox(octree.root()[7][7]) == CGAL::Bbox_3(0.75, 0.75, 0.75, 1.5, 1.5, 1.5)); + assert(octree.bbox(octree.node(7, 0)) == + Octree::Bbox(0, 0, 0, 0.5, 0.5, 0.5)); + assert(octree.bbox(octree.node(7, 1)) == + Octree::Bbox(0.5, 0, 0, 1, 0.5, 0.5)); + assert(octree.bbox(octree.node(7, 2)) == + Octree::Bbox(0, 0.5, 0, 0.5, 1, 0.5)); + assert(octree.bbox(octree.node(7, 3)) == + Octree::Bbox(0.5, 0.5, 0, 1, 1, 0.5)); + assert(octree.bbox(octree.node(7, 4)) == + Octree::Bbox(0, 0, 0.5, 0.5, 0.5, 1)); + assert(octree.bbox(octree.node(7, 5)) == + Octree::Bbox(0.5, 0, 0.5, 1, 0.5, 1)); + assert(octree.bbox(octree.node(7, 6)) == + Octree::Bbox(0, 0.5, 0.5, 0.5, 1, 1)); + assert(octree.bbox(octree.node(7, 7)) == + Octree::Bbox(0.5, 0.5, 0.5, 1, 1, 1)); + + // All child nodes should share a vertex + auto center_of_last_child = octree.bbox(octree.node(7, 7)).vertex(0); + assert(octree.bbox(octree.node(7, 0)).vertex(7) == center_of_last_child); + assert(octree.bbox(octree.node(7, 1)).vertex(4) == center_of_last_child); + assert(octree.bbox(octree.node(7, 2)).vertex(6) == center_of_last_child); + assert(octree.bbox(octree.node(7, 3)).vertex(5) == center_of_last_child); + assert(octree.bbox(octree.node(7, 4)).vertex(2) == center_of_last_child); + assert(octree.bbox(octree.node(7, 5)).vertex(3) == center_of_last_child); + assert(octree.bbox(octree.node(7, 6)).vertex(1) == center_of_last_child); + + // Nodes of different sizes should also share vertices + assert(octree.bbox(octree.node(7, 0)).vertex(0) == octree.bbox(octree.node(0, 7)).vertex(7)); } int main(void) { diff --git a/Orthtree/test/Orthtree/test_octree_copy_move_constructors.cpp b/Orthtree/test/Orthtree/test_octree_copy_move_constructors.cpp index ae6272fe149e..e802572ccff4 100644 --- a/Orthtree/test/Orthtree/test_octree_copy_move_constructors.cpp +++ b/Orthtree/test/Orthtree/test_octree_copy_move_constructors.cpp @@ -1,51 +1,65 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include -#include +#include +#include + #include +#include +#include +#include #include -#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; +using Octree_without_data = CGAL::Orthtree>; -int main(void) +template +int test(Tree &tree) { - std::size_t nb_pts = 100; - Point_set points; - CGAL::Random_points_in_cube_3 generator; - points.reserve(nb_pts); - for (std::size_t i = 0; i < nb_pts; ++i) - points.insert(*(generator++)); + assert (tree.is_leaf(tree.root())); // tree is not refined yet - Octree base (points, points.point_map()); - assert (base.root().is_leaf()); // base is not refined yet + Tree copy1 (tree); + assert (copy1.is_leaf(copy1.root())); // copy1 is thus not refined either + assert (tree == copy1); // tree should be equal to copy1 - Octree copy1 (base); - assert (copy1.root().is_leaf()); // copy1 is thus not refined either - assert (base == copy1); // base should be equal to copy1 + tree.refine(CGAL::Orthtrees::Maximum_depth(5)); + assert (!tree.is_leaf(tree.root())); // tree is now refined + assert (copy1.is_leaf(copy1.root())); // copy1 should be unaffected and still unrefined + assert (tree != copy1); // tree should be different from copy1 - base.refine(); - assert (!base.root().is_leaf()); // base is now refined - assert (copy1.root().is_leaf()); // copy1 should be unaffected and still unrefined - assert (base != copy1); // base should be different from copy1 + Tree copy2 (tree); + assert (!copy2.is_leaf(copy2.root())); // copy2 should be refined + assert (tree == copy2); // tree should be equal to copy2 - Octree copy2 (base); - assert (!copy2.root().is_leaf()); // copy2 should be refined - assert (base == copy2); // base should be equal to copy2 + Tree move (std::move(tree)); + assert (!move.is_leaf(move.root())); // move should be refined - Octree move (std::move(base)); - assert (!move.root().is_leaf()); // move should be refined - assert (base.root().is_leaf()); // base should be back to init state (unrefined) - assert (copy1.root().is_leaf()); // copy1 still unaffected and still unrefined - assert (!copy2.root().is_leaf()); // copy2 unaffected by move and still refined + assert (tree.is_leaf(tree.root())); // tree should be back to init state (unrefined) + assert (copy1.is_leaf(copy1.root())); // copy1 still unaffected and still unrefined + assert (!copy2.is_leaf(copy2.root())); // copy2 unaffected by move and still refined assert (move == copy2); // move should be equal to copy2 return EXIT_SUCCESS; } + +int main() +{ + std::size_t nb_pts = 100; + Point_set points; + CGAL::Random_points_in_cube_3 generator; + points.reserve(nb_pts); + for (std::size_t i = 0; i < nb_pts; ++i) + points.insert(*(generator++)); + + Octree base(points, points.point_map()); + test(base); + + Octree_without_data base2({}); + test(base2); +} diff --git a/Orthtree/test/Orthtree/test_octree_custom_properties.cpp b/Orthtree/test/Orthtree/test_octree_custom_properties.cpp new file mode 100644 index 000000000000..2c4ab94918b1 --- /dev/null +++ b/Orthtree/test/Orthtree/test_octree_custom_properties.cpp @@ -0,0 +1,72 @@ +#define CGAL_TRACE_STREAM std::cerr + +#include + +#include +#include +#include + +#include +#include + +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; + +int main(void) { + + std::size_t nb_pts = 100; + Point_set points; + CGAL::Random_points_in_cube_3 generator; + points.reserve(nb_pts); + for (std::size_t i = 0; i < nb_pts; ++i) + points.insert(*(generator++)); + + Octree tree(points, points.point_map()); + + // Testing built in node properties + typename Octree::Property_map data_prop = *tree.property("contents"); + CGAL_USE(data_prop); + + // list of properties + assert(tree.properties().size() == 5); + auto prop1 = tree.add_property("test", int(5)); + assert(prop1.second); + // One property added + assert(tree.properties().size() == 6); + // Default value should be respected + assert(prop1.first[tree.root()] == 5); + // Changes to individual nodes should be respected + prop1.first[tree.root()] = 0; + assert(prop1.first[tree.root()] == 0); + + auto prop2 = tree.add_property("test", int(0)); + assert(!prop2.second); + assert(tree.properties().size() == 6); + + auto prop3 = tree.add_property("test2", std::string()); + assert(prop3.second); + assert(tree.properties().size() == 7); + + auto prop4 = tree.property("test"); + assert(prop4.has_value()); + + auto prop5 = tree.property("test"); + assert(!prop5.has_value()); + + // Expanding the tree; new nodes should be assigned the default value + tree.refine(10, 1); + for (auto n : tree.traverse>()) { + // Everything but the root will have the default value + if (!tree.is_root(n)) assert(prop1.first[n] == 5); + } + // The root should have preserved its custom value + assert(prop1.first[tree.root()] == 0); + + tree.remove_property(prop1.first); + assert(tree.properties().size() == 6); + + return EXIT_SUCCESS; +} diff --git a/Orthtree/test/Orthtree/test_octree_equality.cpp b/Orthtree/test/Orthtree/test_octree_equality.cpp index 1eb88263774e..ce6d3c95dcb6 100644 --- a/Orthtree/test/Orthtree/test_octree_equality.cpp +++ b/Orthtree/test/Orthtree/test_octree_equality.cpp @@ -1,16 +1,15 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include -#include #include +#include +#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; void test_identical_trees() { diff --git a/Orthtree/test/Orthtree/test_octree_grade.cpp b/Orthtree/test/Orthtree/test_octree_grade.cpp index a0a574bd9eff..19acd787d9a0 100644 --- a/Orthtree/test/Orthtree/test_octree_grade.cpp +++ b/Orthtree/test/Orthtree/test_octree_grade.cpp @@ -1,36 +1,36 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include #include -#include -#include -#include +#include #include +#include +#include +#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree Octree; -typedef CGAL::Orthtrees::Leaves_traversal Leaves_traversal; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; +using Leaves_traversal = CGAL::Orthtrees::Leaves_traversal; -std::size_t count_jumps(Octree &octree) { +std::size_t count_jumps(Octree& octree) { std::size_t jumps = 0; - for (auto &node : octree.traverse(Leaves_traversal())) { + for (auto node: octree.traverse()) { for (int direction = 0; direction < 6; ++direction) { - auto adjacent_node = node.adjacent_node(direction); + auto adjacent_node = octree.adjacent_node(node, direction); - if (adjacent_node.is_null()) + if (!adjacent_node) continue; - if ((node.depth() - adjacent_node.depth()) > 1) + if ((octree.depth(node) - octree.depth(*adjacent_node)) > 1) jumps++; } } diff --git a/Orthtree/test/Orthtree/test_octree_intersecting.cpp b/Orthtree/test/Orthtree/test_octree_intersecting.cpp index 59ce353fd28a..7576ad168ec4 100644 --- a/Orthtree/test/Orthtree/test_octree_intersecting.cpp +++ b/Orthtree/test/Orthtree/test_octree_intersecting.cpp @@ -1,17 +1,17 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include #include -#include -#include #include +#include +#include +#include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef std::vector Point_vector; -typedef CGAL::Octree Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_vector = std::vector; +using Octree = CGAL::Octree; int main(void) { @@ -34,7 +34,7 @@ int main(void) { points.emplace_back(-0.9, -1, -1); // Create an octree from the vector - Octree octree(points); + Octree octree(Octree::Traits{points}); // Build the octree octree.refine(10, 2); @@ -46,7 +46,7 @@ int main(void) { auto query = Point{1, 1, 1}; // Get a list of nodes intersected - std::vector nodes{}; + std::vector nodes{}; octree.intersected_nodes(query, std::back_inserter(nodes)); // A point should only intersect one node @@ -63,15 +63,15 @@ int main(void) { auto query = Kernel::Sphere_3(Point{1, 0.5, 1}, 1.0); // Get a list of nodes intersected - std::vector nodes{}; + std::vector nodes{}; octree.intersected_nodes(query, std::back_inserter(nodes)); // Check the results assert(4 == nodes.size()); - assert(octree[Octree::Traits::RIGHT_TOP_BACK] == nodes[0]); - assert(octree[Octree::Traits::RIGHT_BOTTOM_FRONT] == nodes[1]); - assert(octree[Octree::Traits::LEFT_TOP_FRONT] == nodes[2]); - assert(octree[Octree::Traits::RIGHT_TOP_FRONT] == nodes[3]); + assert(octree.node(Octree::Traits::RIGHT_TOP_BACK) == nodes[0]); + assert(octree.node(Octree::Traits::RIGHT_BOTTOM_FRONT) == nodes[1]); + assert(octree.node(Octree::Traits::LEFT_TOP_FRONT) == nodes[2]); + assert(octree.node(Octree::Traits::RIGHT_TOP_FRONT) == nodes[3]); } // Intersection with a ray @@ -81,19 +81,22 @@ int main(void) { auto query = Kernel::Ray_3(Point{1, 1, 1}, Point{0, 0, 0}); // Get a list of nodes intersected - std::vector nodes{}; + std::vector nodes{}; octree.intersected_nodes(query, std::back_inserter(nodes)); // Check the results assert(8 == nodes.size()); - assert(octree[Octree::Traits::LEFT_BOTTOM_BACK] == nodes[0]); - assert(octree[Octree::Traits::RIGHT_BOTTOM_BACK][Octree::Traits::LEFT_TOP_FRONT] == nodes[1]); - assert(octree[Octree::Traits::LEFT_TOP_BACK] == nodes[2]); - assert(octree[Octree::Traits::RIGHT_TOP_BACK] == nodes[3]); - assert(octree[Octree::Traits::LEFT_BOTTOM_FRONT] == nodes[4]); - assert(octree[Octree::Traits::RIGHT_BOTTOM_FRONT] == nodes[5]); - assert(octree[Octree::Traits::LEFT_TOP_FRONT] == nodes[6]); - assert(octree[Octree::Traits::RIGHT_TOP_FRONT] == nodes[7]); + assert(octree.node(Octree::Traits::LEFT_BOTTOM_BACK) == nodes[0]); + assert( + octree.node(Octree::Traits::RIGHT_BOTTOM_BACK, Octree::Traits::LEFT_TOP_FRONT) + == nodes[1] + ); + assert(octree.node(Octree::Traits::LEFT_TOP_BACK) == nodes[2]); + assert(octree.node(Octree::Traits::RIGHT_TOP_BACK) == nodes[3]); + assert(octree.node(Octree::Traits::LEFT_BOTTOM_FRONT) == nodes[4]); + assert(octree.node(Octree::Traits::RIGHT_BOTTOM_FRONT) == nodes[5]); + assert(octree.node(Octree::Traits::LEFT_TOP_FRONT) == nodes[6]); + assert(octree.node(Octree::Traits::RIGHT_TOP_FRONT) == nodes[7]); } return EXIT_SUCCESS; diff --git a/Orthtree/test/Orthtree/test_octree_kernels.cpp b/Orthtree/test/Orthtree/test_octree_kernels.cpp index 58d49eb4fccd..5e477d4a9ebd 100644 --- a/Orthtree/test/Orthtree/test_octree_kernels.cpp +++ b/Orthtree/test/Orthtree/test_octree_kernels.cpp @@ -1,16 +1,16 @@ #include +#include +#include #include #include #include -#include -#include template void test() { - typedef typename Kernel::Point_3 Point; - typedef CGAL::Point_set_3 Point_set; - typedef CGAL::Octree Octree; + using Point = typename Kernel::Point_3; + using Point_set = CGAL::Point_set_3; + using Octree = CGAL::Octree; Point_set points; CGAL::Random_points_in_cube_3 generator; diff --git a/Orthtree/test/Orthtree/test_octree_locate.cpp b/Orthtree/test/Orthtree/test_octree_locate.cpp index d601e3614660..8bbc8932d9ea 100644 --- a/Orthtree/test/Orthtree/test_octree_locate.cpp +++ b/Orthtree/test/Orthtree/test_octree_locate.cpp @@ -1,19 +1,18 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include -#include #include +#include +#include #include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; void test_1_point() { @@ -52,24 +51,24 @@ void test_8_points() { octree.refine(10, 1); // Existing points should end up in the same place - assert(octree.root()[0] == octree.locate({-1, -1, -1})); - assert(octree.root()[1] == octree.locate({1, -1, -1})); - assert(octree.root()[2] == octree.locate({-1, 1, -1})); - assert(octree.root()[3] == octree.locate({1, 1, -1})); - assert(octree.root()[4] == octree.locate({-1, -1, 1})); - assert(octree.root()[5] == octree.locate({1, -1, 1})); - assert(octree.root()[6] == octree.locate({-1, 1, 1})); - assert(octree.root()[7] == octree.locate({1, 1, 1})); + assert(octree.node(0) == octree.locate({-1, -1, -1})); + assert(octree.node(1) == octree.locate({1, -1, -1})); + assert(octree.node(2) == octree.locate({-1, 1, -1})); + assert(octree.node(3) == octree.locate({1, 1, -1})); + assert(octree.node(4) == octree.locate({-1, -1, 1})); + assert(octree.node(5) == octree.locate({1, -1, 1})); + assert(octree.node(6) == octree.locate({-1, 1, 1})); + assert(octree.node(7) == octree.locate({1, 1, 1})); // Points adjacent to the existing points should also end up in the same place - assert(octree.root()[0] == octree.locate({-1.1, -1.1, -1.1})); - assert(octree.root()[1] == octree.locate({1.1, -1.1, -1.1})); - assert(octree.root()[2] == octree.locate({-1.1, 1.1, -1.1})); - assert(octree.root()[3] == octree.locate({1.1, 1.1, -1.1})); - assert(octree.root()[4] == octree.locate({-1.1, -1.1, 1.1})); - assert(octree.root()[5] == octree.locate({1.1, -1.1, 1.1})); - assert(octree.root()[6] == octree.locate({-1.1, 1.1, 1.1})); - assert(octree.root()[7] == octree.locate({1.1, 1.1, 1.1})); + assert(octree.node(0) == octree.locate({-0.99, -0.99, -0.99})); + assert(octree.node(1) == octree.locate({0.99, -0.99, -0.99})); + assert(octree.node(2) == octree.locate({-0.99, 0.99, -0.99})); + assert(octree.node(3) == octree.locate({0.99, 0.99, -0.99})); + assert(octree.node(4) == octree.locate({-0.99, -0.99, 0.99})); + assert(octree.node(5) == octree.locate({0.99, -0.99, 0.99})); + assert(octree.node(6) == octree.locate({-0.99, 0.99, 0.99})); + assert(octree.node(7) == octree.locate({0.99, 0.99, 0.99})); } @@ -93,24 +92,24 @@ void test_10_points() { octree.refine(10, 1); // Existing points should end up in the same place - assert(octree.root()[0] == octree.locate({-1, -1, -1})); - assert(octree.root()[1] == octree.locate({1, -1, -1})); - assert(octree.root()[2] == octree.locate({-1, 1, -1})); - assert(octree.root()[3][3][3] == octree.locate({1, 1, -1})); - assert(octree.root()[4][4][4] == octree.locate({-1, -1, 1})); - assert(octree.root()[5] == octree.locate({1, -1, 1})); - assert(octree.root()[6] == octree.locate({-1, 1, 1})); - assert(octree.root()[7] == octree.locate({1, 1, 1})); + assert(octree.node(0) == octree.locate({-1, -1, -1})); + assert(octree.node(1) == octree.locate({1, -1, -1})); + assert(octree.node(2) == octree.locate({-1, 1, -1})); + assert(octree.node(3, 3, 3, 3, 3) == octree.locate({1, 1, -1})); + assert(octree.node(4, 4, 4) == octree.locate({-1, -1, 1})); + assert(octree.node(5) == octree.locate({1, -1, 1})); + assert(octree.node(6) == octree.locate({-1, 1, 1})); + assert(octree.node(7) == octree.locate({1, 1, 1})); // Points adjacent to the existing points might end up in different places - assert(octree.root()[0] == octree.locate({-1.1, -1.1, -1.1})); - assert(octree.root()[1] == octree.locate({1.1, -1.1, -1.1})); - assert(octree.root()[2] == octree.locate({-1.1, 1.1, -1.1})); - assert(octree.root()[3][3][3] == octree.locate({1.1, 1.1, -1.1})); - assert(octree.root()[4][4][4] == octree.locate({-1.1, -1.1, 1.1})); - assert(octree.root()[5] == octree.locate({1.1, -1.1, 1.1})); - assert(octree.root()[6] == octree.locate({-1.1, 1.1, 1.1})); - assert(octree.root()[7] == octree.locate({1.1, 1.1, 1.1})); + assert(octree.node(0) == octree.locate({-0.99, -0.99, -0.99})); + assert(octree.node(1) == octree.locate({0.99, -0.99, -0.99})); + assert(octree.node(2) == octree.locate({-0.99, 0.99, -0.99})); + assert(octree.node(3, 3, 3, 3, 3) == octree.locate({0.99, 0.99, -0.99})); + assert(octree.node(4, 4, 4) == octree.locate({-0.99, -0.99, 0.99})); + assert(octree.node(5) == octree.locate({0.99, -0.99, 0.99})); + assert(octree.node(6) == octree.locate({-0.99, 0.99, 0.99})); + assert(octree.node(7) == octree.locate({0.99, 0.99, 0.99})); } diff --git a/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp b/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp index 8e87eaab152d..2f37a75236cb 100644 --- a/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp +++ b/Orthtree/test/Orthtree/test_octree_nearest_neighbor.cpp @@ -1,30 +1,28 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include -#include #include #include #include #include #include +#include +#include #include #include using namespace std::chrono; -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::FT FT; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree -Octree; - -typedef CGAL::Search_traits_3 Kd_tree_traits; -typedef CGAL::Orthogonal_k_neighbor_search Kd_tree_search; -typedef Kd_tree_search::Tree Kd_tree; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using FT = Kernel::FT; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; +using Kd_tree_traits = CGAL::Search_traits_3; +using Kd_tree_search = CGAL::Orthogonal_k_neighbor_search; +using Kd_tree = Kd_tree_search::Tree; void naive_vs_octree(std::size_t dataset_size) { @@ -47,7 +45,7 @@ void naive_vs_octree(std::size_t dataset_size) { { FT distance_nearest = (std::numeric_limits::max)(); - for (auto &p : points.points()) { + for (auto& p: points.points()) { FT distance_current = CGAL::squared_distance(p, random_point); if (distance_current < distance_nearest) { @@ -72,10 +70,9 @@ void naive_vs_octree(std::size_t dataset_size) { octree.refine(10, 20); auto octree_start_time = high_resolution_clock::now(); { - // TODO: Write a nearest-neighbor implementation and use it here - std::vector k_neighbors; - octree.nearest_neighbors(random_point, 1, std::back_inserter(k_neighbors)); - octree_nearest = *k_neighbors.begin(); + std::vector k_neighbors; + octree.nearest_k_neighbors(random_point, 1, std::back_inserter(k_neighbors)); + octree_nearest = get(points.point_map(), *k_neighbors.begin()); } duration octree_elapsed_time = high_resolution_clock::now() - octree_start_time; @@ -109,9 +106,9 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) { Kd_tree kd_tree(points.points().begin(), points.points().end()); kd_tree.build(); auto kd_tree_start_time = high_resolution_clock::now(); - Kd_tree_search search(kd_tree, random_point, (unsigned int)(K)); + Kd_tree_search search(kd_tree, random_point, (unsigned int) (K)); duration kd_tree_elapsed_time = high_resolution_clock::now() - kd_tree_start_time; - for (auto p : search) + for (auto p: search) kd_tree_nearest_neighbors.push_back(p.first); std::cout << "Kd_tree --> " @@ -120,11 +117,11 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) { << std::endl; // Do the same using the octree - std::vector octree_nearest_neighbors; + std::vector octree_nearest_neighbors; Octree octree(points, points.point_map()); octree.refine(10, 20); auto octree_start_time = high_resolution_clock::now(); - octree.nearest_neighbors(random_point, K, std::back_inserter(octree_nearest_neighbors)); + octree.nearest_k_neighbors(random_point, K, std::back_inserter(octree_nearest_neighbors)); duration octree_elapsed_time = high_resolution_clock::now() - octree_start_time; std::cout << "Octree --> " @@ -137,12 +134,13 @@ void kdtree_vs_octree(std::size_t dataset_size, std::size_t K) { // Check that they produce the same answer for (std::size_t j = 0; j < K; ++j) - assert(octree_nearest_neighbors[j] == kd_tree_nearest_neighbors[j]); + assert(get(points.point_map(), octree_nearest_neighbors[j]) == kd_tree_nearest_neighbors[j]); } int main(void) { + naive_vs_octree(21); naive_vs_octree(500); naive_vs_octree(1000); naive_vs_octree(10000); diff --git a/Orthtree/test/Orthtree/test_octree_refine.cpp b/Orthtree/test/Orthtree/test_octree_refine.cpp index ff5177821790..ff75dd18edc7 100644 --- a/Orthtree/test/Orthtree/test_octree_refine.cpp +++ b/Orthtree/test/Orthtree/test_octree_refine.cpp @@ -1,18 +1,34 @@ #define CGAL_TRACE_STREAM std::cerr -#include #include -#include #include +#include +#include #include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree Octree; -typedef Octree::Node Node; +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; + +class Split_nth_child_of_root { + std::size_t m_n; +public: + + Split_nth_child_of_root(std::size_t n) : m_n(n) {} + + template + bool operator()(const Node& node) const { + return (node.depth() == 1 && node.local_coordinates().to_ulong() == m_n); + } + + template + bool operator()(Node_index i, const Tree &tree) const { + return (tree.depth(i) == 1 && tree.local_coordinates(i).to_ulong() == m_n); + } +}; void test_1_point() { @@ -24,13 +40,9 @@ void test_1_point() { Octree octree(points, points.point_map()); octree.refine(10, 1); - // Check that the topology matches - Node single_node = CGAL::Orthtrees::Node_access::create_node(Node(), 0); - CGAL::Orthtrees::Node_access::points(single_node) - = CGAL::Orthtrees::Node_access::points(octree.root()); - assert(Node::is_topology_equal(single_node, octree.root())); + // Check that the root node was never split + assert(octree.is_leaf(octree.root())); assert(0 == octree.depth()); - CGAL::Orthtrees::Node_access::free(single_node); } void test_2_points() { @@ -45,12 +57,10 @@ void test_2_points() { octree.refine(10, 1); // The octree should have been split once - Node other = CGAL::Orthtrees::Node_access::create_node(Node(), 0); - CGAL::Orthtrees::Node_access::split(other); - assert(Node::is_topology_equal(other, octree.root())); + Octree other(points, points.point_map()); + other.split(other.root()); + assert(Octree::is_topology_equal(other, octree)); assert(1 == octree.depth()); - CGAL::Orthtrees::Node_access::free(other); - } void test_4_points() { @@ -65,14 +75,18 @@ void test_4_points() { Octree octree(points, points.point_map()); octree.refine(10, 1); - // The octree should have been split once on the first level, and twice on the second - Node other = CGAL::Orthtrees::Node_access::create_node(Node(), 0); - CGAL::Orthtrees::Node_access::split(other); - CGAL::Orthtrees::Node_access::split(other[3]); - CGAL::Orthtrees::Node_access::split(other[7]); - assert(Node::is_topology_equal(other, octree.root())); + Octree other(points, points.point_map()); + other.split(other.root()); + other.split(other.node(3)); + other.split(other.node(7)); + assert(Octree::is_topology_equal(other, octree)); assert(2 == octree.depth()); - CGAL::Orthtrees::Node_access::free(other); + + // Applying another splitting criterion shouldn't reset the tree. + octree.refine(Split_nth_child_of_root(2)); + other.split(other.node(2)); + assert(Octree::is_topology_equal(other, octree)); + } int main(void) { diff --git a/Orthtree/test/Orthtree/test_octree_traverse.cpp b/Orthtree/test/Orthtree/test_octree_traverse.cpp index 57853df8dfaf..e7ec64fd3577 100644 --- a/Orthtree/test/Orthtree/test_octree_traverse.cpp +++ b/Orthtree/test/Orthtree/test_octree_traverse.cpp @@ -2,17 +2,18 @@ #include #include - -#include #include +#include #include -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef CGAL::Point_set_3 Point_set; -typedef CGAL::Octree Octree; -typedef CGAL::Orthtrees::Preorder_traversal Preorder_traversal; + +using Kernel = CGAL::Simple_cartesian; +using Point = Kernel::Point_3; +using Point_set = CGAL::Point_set_3; +using Octree = CGAL::Octree; +using Preorder_traversal = CGAL::Orthtrees::Preorder_traversal; +using Level_traversal = CGAL::Orthtrees::Level_traversal; bool test_preorder_1_node() { @@ -53,7 +54,31 @@ bool test_preorder_9_nodes() { assert(*iter == octree.root()); for (int i = 0; i < 8; ++i) { iter++; - assert((*iter == octree.root()[i])); + assert(*iter == octree.node(i)); + } + + return true; +} + +bool test_level_9_nodes() { + + // Define the dataset + Point_set points; + points.insert({-1, -1, -1}); + points.insert({1, -1, -1}); + + // Create the octree + Octree octree(points, points.point_map()); + octree.refine(10, 1); + + // Create the range + auto nodes = octree.traverse(static_cast(1)); + + // Check each item in the range + auto iter = nodes.begin(); + for (int i = 0; i < 8; ++i) { + assert(*iter == octree.node(i)); + iter++; } return true; @@ -71,6 +96,7 @@ bool test_preorder_25_nodes() { // Create the octree Octree octree(points, points.point_map()); octree.refine(10, 1); + std::cout << octree << std::endl; // Create the range auto nodes = octree.traverse(); @@ -79,28 +105,28 @@ bool test_preorder_25_nodes() { auto iter = nodes.begin(); assert(*iter == octree.root()); iter++; - assert((*iter == octree.root()[0])); + assert(*iter == octree.node(0)); iter++; - assert((*iter == octree.root()[1])); + assert(*iter == octree.node(1)); iter++; - assert((*iter == octree.root()[2])); + assert((*iter == octree.node(2))); iter++; - assert((*iter == octree.root()[3])); + assert(*iter == octree.node(3)); for (int i = 0; i < 8; ++i) { iter++; - assert((*iter == octree.root()[3][i])); + assert(*iter == octree.node(3, i)); } iter++; - assert((*iter == octree.root()[4])); + assert((*iter == octree.node(4))); iter++; - assert((*iter == octree.root()[5])); + assert((*iter == octree.node(5))); iter++; - assert((*iter == octree.root()[6])); + assert((*iter == octree.node(6))); iter++; - assert((*iter == octree.root()[7])); + assert((*iter == octree.node(7))); for (int i = 0; i < 8; ++i) { iter++; - assert((*iter == octree.root()[7][i])); + assert(*iter == octree.node(7, i)); } return true; @@ -110,6 +136,7 @@ int main(void) { test_preorder_1_node(); test_preorder_9_nodes(); + test_level_9_nodes(); test_preorder_25_nodes(); return 0; diff --git a/Property_map/include/CGAL/Property_container.h b/Property_map/include/CGAL/Property_container.h new file mode 100644 index 000000000000..739f652b3184 --- /dev/null +++ b/Property_map/include/CGAL/Property_container.h @@ -0,0 +1,630 @@ +// Copyright (c) 2023 INRIA +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Jackson Campolattaro + +#ifndef CGAL_PROPERTY_CONTAINTER_H +#define CGAL_PROPERTY_CONTAINTER_H + +#include + +#include +#include +#include + +#include + +#ifndef DOXYGEN_RUNNING + +namespace CGAL::Properties::Experimental { + +template +class Property_array_base { +public: + + Property_array_base() = default; + + Property_array_base(const Property_array_base& rhs) = delete; + + virtual ~Property_array_base() = default; + + // Declare virtual functions here, for things which need to be done within the Property container + // todo: these should mostly be private, and made available using friend + + virtual std::shared_ptr> empty_clone(const std::vector& active_indices) = 0; + + virtual std::shared_ptr> clone(const std::vector& active_indices) = 0; + + virtual void copy(const Property_array_base& other) = 0; + + // desactived as MSVC 2017 as an issue with that but it is not currently used. +#if 0 + virtual void move(Property_array_base&& other) = 0; +#endif + + virtual void append(const Property_array_base& other) = 0; + + virtual void reserve(std::size_t n) = 0; + + virtual void shrink_to_fit() = 0; + + virtual void swap(Index a, Index b) = 0; + + virtual void reset(Index i) = 0; + + virtual const std::type_info& type() const = 0; + + virtual void transfer_from(const Property_array_base& other_base, Index other_index, Index this_index) = 0; + +}; + +/*! + * \brief Indexed storage for arbitrary types + * + * todo: make this effectively private, prioritize the use of Property_array_handle + * + * @tparam T + */ +template +class Property_array : public Property_array_base { + + std::vector m_data; + const std::vector& m_active_indices; + T m_default_value; + +public: + + using value_type = T; + using reference = typename std::vector::reference; + using const_reference = typename std::vector::const_reference; + using iterator = typename std::vector::iterator; + using const_iterator = typename std::vector::const_iterator; + + Property_array(const std::vector& active_indices, const T& default_value) : + m_data(), m_active_indices(active_indices), m_default_value(default_value) { + + m_data.reserve(active_indices.capacity()); + m_data.resize(active_indices.size(), m_default_value); + } + + virtual std::shared_ptr> empty_clone(const std::vector& active_indices) override { + return std::make_shared>(active_indices, m_default_value); + } + + virtual std::shared_ptr> clone(const std::vector& active_indices) override { + auto new_array = std::make_shared>(active_indices, m_default_value); + new_array->m_data = m_data; + return new_array; + } + + virtual void copy(const Property_array_base& other_base) override { + auto& other = dynamic_cast&>(other_base); + m_data = other.m_data; + CGAL_precondition(m_active_indices.size() == m_data.size()); + } + +// deactived as MSVC 2017 as an issue with that but it is not currently used. +#if 0 + virtual void move(Property_array_base&& other_base) override { + auto&& other = static_cast&&>(other_base); + m_data = std::move(other.m_data); + CGAL_precondition(m_active_indices.size() == m_data.size()); + } +#endif + + virtual void append(const Property_array_base& other_base) override { + auto& other = dynamic_cast&>(other_base); + CGAL_precondition(m_data.size() + other.m_data.size() == m_active_indices.size()); + m_data.insert(m_data.end(), other.m_data.begin(), other.m_data.end()); + } + + virtual void reserve(std::size_t n) override { + CGAL_precondition(m_active_indices.size() == n); + m_data.resize(n, m_default_value); + }; + + virtual void shrink_to_fit() override { + m_data.shrink_to_fit(); + } + + virtual void swap(Index a, Index b) override { + // todo: maybe cast to index, instead of casting index to size? + CGAL_precondition(std::size_t(a) < m_data.size() && std::size_t(b) < m_data.size()); + std::iter_swap(m_data.begin() + a, m_data.begin() + b); + }; + + virtual void reset(Index i) override { + CGAL_precondition(std::size_t(i) < m_data.size()); + m_data[std::size_t(i)] = m_default_value; + }; + + virtual const std::type_info& type() const override { return typeid(T); }; + + virtual void transfer_from(const Property_array_base& other_base, + Index other_index, Index this_index) override { + + CGAL_precondition(other_base.type() == type()); + auto& other = dynamic_cast&>(other_base); + CGAL_precondition(std::size_t(other_index) < other.capacity() && std::size_t(this_index) < capacity()); + m_data[this_index] = other.m_data[other_index]; + } + +public: + + // todo: there's not really a good reason to use these, maybe they should be removed + + [[nodiscard]] std::size_t size() const { return std::count(m_active_indices.begin(), m_active_indices.end(), true); } + + [[nodiscard]] std::size_t capacity() const { return m_data.size(); } + + const_reference operator[](Index i) const { + CGAL_precondition(std::size_t(i) < m_data.size()); + return m_data[std::size_t(i)]; + } + + reference operator[](Index i) { + CGAL_precondition(std::size_t(i) < m_data.size()); + return m_data[std::size_t(i)]; + } + + iterator begin() { return m_data.begin(); } + + iterator end() { return m_data.end(); } + + const_iterator begin() const { return m_data.begin(); } + + const_iterator end() const { return m_data.end(); } + +public: + + bool operator==(const Property_array& other) const { + return &other == this; + } + + bool operator!=(const Property_array& other) const { return !operator==(other); } + +}; + +// todo: property maps/array handles should go in their own file + +// todo: add const/read-only handle +template +class Property_array_handle { + + std::reference_wrapper> m_array; + +public: + + // Necessary for use as a boost::property_type + using key_type = Index; + using value_type = T; + using reference = typename std::vector::reference; + using const_reference = typename std::vector::const_reference; + using category = boost::lvalue_property_map_tag; + + using iterator = typename std::vector::iterator; + using const_iterator = typename std::vector::const_iterator; + + Property_array_handle(Property_array& array) : m_array(array) {} + + [[nodiscard]] std::size_t size() const { return m_array.get().size(); } + + [[nodiscard]] std::size_t capacity() const { return m_array.get().capacity(); } + + Property_array& array() const { return m_array.get(); } + + // todo: This might not be needed, if the other operator[] is made const + const_reference operator[](Index i) const { return m_array.get()[i]; } + + reference operator[](Index i) { return m_array.get()[i]; } + + // todo: maybe these can be const, in an lvalue property map? + iterator begin() { return m_array.get().begin(); } + + iterator end() { return m_array.get().end(); } + + const_iterator begin() const { return m_array.get().begin(); } + + const_iterator end() const { return m_array.get().end(); } + + bool operator==(const Property_array_handle& other) const { return other.m_array.get() == m_array.get(); } + + bool operator!=(const Property_array_handle& other) const { return !operator==(other); } + + inline friend reference get(Property_array_handle p, const Index& i) { return p[i]; } + + inline friend void put(Property_array_handle p, const Index& i, const T& v) { p[i] = v; } + +}; + +template +class Property_container { + + std::multimap>> m_properties; + std::vector m_active_indices{}; + +public: + + template + using Array = Property_array; + + Property_container() = default; + + Property_container(const Property_container& other) { + m_active_indices = other.m_active_indices; + + for (auto [name, array] : other.m_properties) { + // todo: this could probably be made faster using emplace_hint + m_properties.emplace( + name, + array->clone(m_active_indices) + ); + } + } + + Property_container(Property_container&& other) { *this = std::move(other); } + + // This is not exactly an assignment as existing unique properties are kept. + Property_container& operator=(const Property_container& other) { + m_active_indices = other.m_active_indices; + + for (auto [name, array] : other.m_properties) { + // search if property already exists + auto range = m_properties.equal_range(name); + auto it = range.first; + for (; it != range.second; it++) { + if (typeid(*array) == typeid((*it->second))) + break; + } + + if (it != range.second) + it->second->copy(*array); + else + m_properties.emplace(name, array->clone(m_active_indices)); + } + + return *this; + } + + // This is not exactly an assignment as existing unique properties are kept. + Property_container& operator=(Property_container&& other) { + m_active_indices = std::move(other.m_active_indices); + + for (auto [name, array] : other.m_properties) { + // search if property already exists + auto range = m_properties.equal_range(name); + auto it = range.first; + for (; it != range.second; it++) { + if (typeid(*array) == typeid((*it->second))) + break; + } + + if (it != range.second) + it->second->copy(std::move(*array)); + else + m_properties.emplace(name, array->clone(m_active_indices)); + } + + // The moved-from property map should retain all of its properties, but contain 0 elements + other.reserve(0); + return *this; + } + + template + std::pair>, bool> + get_or_add_property(const std::string& name, const T default_value = T()) { + auto range = m_properties.equal_range(name); + for (auto it = range.first; it != range.second; it++) { + Property_array* typed_array_ptr = dynamic_cast*>(it->second.get()); + if (typed_array_ptr != nullptr) + return { {*typed_array_ptr}, false }; + } + + auto it = m_properties.emplace( + name, + std::make_shared>( + m_active_indices, + default_value + ) + ); + + return {{*dynamic_cast*>(it->second.get())}, true}; + } + + template + Property_array& add_property(const std::string& name, const T default_value = T()) { + // todo: I'm not settled on the naming, but it's really convenient to have a function like this + auto [array, created] = get_or_add_property(name, default_value); + CGAL_precondition(created); + return array.get(); + } + +/* + // todo: misleading name, maybe it could be add_same_properties? + void copy_properties(const Property_container& other) { + for (auto [name, other_array]: other.m_properties) { + // If this container doesn't have any property by this name, add it (with the same type as in other) + if (!property_exists(name)) + m_property_arrays.emplace(name, other_array->empty_clone(m_active_indices)); + } + }*/ + + template + const Property_array& get_property(const std::string& name) const { + return *(get_property_if_exists(name)); + } + + template + Property_array& get_property(const std::string& name) { + return *(get_property_if_exists(name)); + } + + template + std::optional>> get_property_if_exists(const std::string& name) const { + auto range = m_properties.equal_range(name); + for (auto it = range.first; it != range.second; it++) { + Property_array* typed_array_ptr = dynamic_cast*>(it->second.get()); + if (typed_array_ptr != nullptr) + return *typed_array_ptr; + } + + return {}; + } + + template + bool property_exists(const std::string& name) const { + auto range = m_properties.equal_range(name); + + for (auto it = range.first; it != range.second; it++) { + Property_array* typed_array_ptr = dynamic_cast*>(it->second.get()); + if (typed_array_ptr != nullptr) + return true; + } + + return false; + } + + /*! + * Removes all properties with the name from the container. + * + * @param name + * @return number of removed properties. + */ + std::size_t remove_properties(const std::string& name) { return m_properties.erase(name); } + + template + bool remove_property(const Property_array& arrayToRemove) { + const Property_array_base* ref = dynamic_cast*>(&arrayToRemove); + for (auto it = m_properties.begin(); it != m_properties.end(); it++) { + auto const& [name, array] = *it; + if (array.get() == ref) { + m_properties.erase(it); + return true; + } + } + return false; + } + + void remove_all_properties_except(const std::vector& preserved_names) { + // todo: if this is used often, it should take a parameter pack instead of a vector + // A fold expression could then be used in place of std::find for better performance + for (auto it = m_properties.begin(); it != m_properties.end();) { + auto const& [name, array] = *it; + if (std::find(preserved_names.begin(), preserved_names.end(), name) == preserved_names.end()) + it = m_properties.erase(it); + else + it++; + } + } + + std::vector properties() const { + std::vector property_names{}; + for (auto const& [name, _]: m_properties) + property_names.emplace_back(name); + return property_names; + } + + std::size_t num_properties() const { return m_properties.size(); } + +/* Deactivated as there may be several Property_maps with different types but the same name. + const std::type_info& property_type(const std::string& name) const { + if (auto it = m_property_arrays.find(name); it != m_property_arrays.end()) + return it->second->type(); + else + return typeid(void); + }*/ + +public: + + void reserve(std::size_t n) { + m_active_indices.resize(n); + for (auto [name, array]: m_properties) + array->reserve(n); + } + + void resize(std::size_t n) { + reserve(n); + std::fill(m_active_indices.begin(), m_active_indices.end(), true); + } + + [[nodiscard]] std::size_t size() const { return std::count(m_active_indices.begin(), m_active_indices.end(), true); } + + [[nodiscard]] std::size_t capacity() const { return m_active_indices.size(); } + + Index emplace_back() { + + // Expand the storage and return the last element + reserve(capacity() + 1); + m_active_indices.back() = true; + auto first_new_index = Index(capacity() - 1); + reset(first_new_index); + return first_new_index; + } + + Index emplace() { + + // If there are empty slots, return the index of one of them and mark it as full + auto first_unused = std::find_if(m_active_indices.begin(), m_active_indices.end(), [](bool used) { return !used; }); + if (first_unused != m_active_indices.end()) { + *first_unused = true; + auto index = Index(std::distance(m_active_indices.begin(), first_unused)); + reset(index); + return index; + } + + return emplace_back(); + } + + Index emplace_group_back(std::size_t n) { + + // Expand the storage and return the start of the new region + reserve(capacity() + n); + for (auto it = m_active_indices.end() - n; it < m_active_indices.end(); ++it) + *it = true; + return Index(capacity() - n); + } + + Index emplace_group(std::size_t n) { + + auto search_start = m_active_indices.begin(); + while (search_start != m_active_indices.end()) { + + // Find the first unused cell + auto unused_begin = std::find_if( + search_start, m_active_indices.end(), + [](bool used) { return !used; } + ); + + auto unused_end = unused_begin; + + // Determine if the group fits + if (std::distance(unused_begin, m_active_indices.end()) >= static_cast::iterator>::difference_type>(n)) + unused_end = std::find_if( + unused_begin, (std::min)(unused_begin + n, m_active_indices.end()), + [](bool used) { return used; } + ); + + // If the discovered range was large enough + if (std::distance(unused_begin, unused_end) >= static_cast::iterator>::difference_type>(n)) { + + // Mark the indices as used, and reset the properties of each of them + // todo: it would be better to provide a function to set a range + for (auto it = unused_begin; it < unused_end; ++it) { + *it = true; + reset(Index(std::distance(m_active_indices.begin(), it))); + } + + // Return the first index of the range + return Index(std::distance(m_active_indices.begin(), unused_begin)); + } + + // If we didn't find a large enough region, continue our search after the end + search_start = unused_end; + } + + // If no empty regions were found, expand the storage + return emplace_group_back(n); + } + + void swap(Index a, Index b) { + for (auto [name, array]: m_properties) + array->swap(a, b); + } + + void reset(Index i) { + for (auto [name, array]: m_properties) + array->reset(i); + } + + void erase(Index i) { + m_active_indices[i] = false; + for (auto [name, array]: m_properties) + array->reset(i); + } + + bool is_erased(Index i) const { + return !m_active_indices[i]; + } + + // todo: I'd prefer to eliminate this, if possible + void mark_active(Index i) { + return m_active_indices[i] = true; + } + + void mark_inactive(Index i) { + return m_active_indices[i] = false; + } + + std::vector active_list() const { + std::vector indices; + for (std::size_t i = 0; i < m_active_indices.size(); ++i) + if (m_active_indices[i]) indices.emplace_back(i); + return indices; + } + + std::vector inactive_list() const { + std::vector indices; + for (std::size_t i = 0; i < m_active_indices.size(); ++i) + if (!m_active_indices[i]) indices.emplace_back(i); + return indices; + } + + void shrink_to_fit() { + for (auto [name, array]: m_properties) + array->shrink_to_fit(); + } + + /*! + * Adds the elements of the other container to this container for each property which is present in this container. + * + * Gaps in both containers are preserved, and all elements of the other container are guaranteed + * to appear after the elements of this container. + * Properties in this container which don't appear in the other container are extended with default values. + * Properties in the other container which don't appear in this one are not included. + * todo: merge() would be useful as well, but could break contiguous regions in the other container + * + * @param other + */ + void append(const Property_container& other) { + + m_active_indices.insert(m_active_indices.end(), other.m_active_indices.begin(), other.m_active_indices.end()); + for (auto [name, array]: m_properties) { + + auto range = other.m_properties.equal_range(name); + auto it = range.first; + for (; it != range.second; it++) { + if (typeid(array.get()) == typeid((it->second.get()))) + break; + } + + if (it != range.second) + array->append(*it->second.get()); + else + array->reserve(m_active_indices.size()); + } + } + +/* + // todo: maybe should be renamed to transfer_from, but I'd rather remove this functionality entirely + void transfer(const Property_container& other, Index other_index, Index this_index) { + CGAL_precondition(other.m_property_arrays.size() == m_property_arrays.size()); + for (auto [name, array]: m_property_arrays) { + auto other_array = other.m_property_arrays.at(name); + array->transfer_from(*other_array, other_index, this_index); + } + }*/ + + // todo: maybe a compress() method? +}; + +} + +#endif // DOXYGEN_RUNNING + +#endif //CGAL_PROPERTY_CONTAINTER_H diff --git a/Property_map/test/Property_map/CMakeLists.txt b/Property_map/test/Property_map/CMakeLists.txt index b0480bfa2098..9f920ea96110 100644 --- a/Property_map/test/Property_map/CMakeLists.txt +++ b/Property_map/test/Property_map/CMakeLists.txt @@ -8,6 +8,7 @@ create_single_source_cgal_program("test_property_map.cpp") create_single_source_cgal_program("dynamic_property_map.cpp") create_single_source_cgal_program("dynamic_properties_test.cpp") create_single_source_cgal_program("kernel_converter_properties_test.cpp") +create_single_source_cgal_program("test_Property_container.cpp") find_package(OpenMesh QUIET) if(OpenMesh_FOUND) diff --git a/Property_map/test/Property_map/test_Property_container.cpp b/Property_map/test/Property_map/test_Property_container.cpp new file mode 100644 index 000000000000..1efe25f4001a --- /dev/null +++ b/Property_map/test/Property_map/test_Property_container.cpp @@ -0,0 +1,261 @@ + +#include +#include + +using namespace CGAL::Properties::Experimental; + +void test_property_creation() { + + Property_container properties; + + // Should return an integer array which didn't previously exist + auto [integers, created] = properties.get_or_add_property("integer", 5); + static_assert(std::is_same_v>>); + assert(created); + assert(properties.num_properties() == 1); + + auto [floats, _] = properties.get_or_add_property("float"); + static_assert(std::is_same_v>>); + assert(properties.num_properties() == 2); + + // get() should retreive the same arrays + assert(integers.get() == properties.get_property("integer")); + assert(floats.get() == properties.get_property("float")); + + // remove() should delete a property array & return if it existed + assert(!properties.remove_properties("not-a-real-property")); + auto removed = properties.remove_property(integers); + assert(removed); + assert(properties.num_properties() == 1); + + // Add a new property + auto [bools, bools_created] = properties.get_or_add_property("bools", false); + static_assert(std::is_same_v>>); + Property_array& b = bools.get(); + CGAL_USE(b); +} + +void test_element_access() { + + Property_container properties; + + auto& integers = properties.add_property("integers", 5); + + // Reserve space for 100 elements + properties.reserve(100); + assert(properties.capacity() == 100); + assert(properties.size() == 0); + + // Newly emplaced elements should go at the front + assert(properties.emplace() == 0); + assert(properties.emplace() == 1); + assert(properties.emplace() == 2); + assert(properties.size() == 3); + + // Make sure that the new elements are equal to the default value + assert(integers[0] == 5); + assert(integers[1] == 5); + assert(integers[2] == 5); + + // Add a new property + auto& floats = properties.add_property("floats", 6.0f); + + // The new property array should already be of the right size + assert(floats.capacity() == 100); + assert(properties.size() == 3); + + // Pre-existing elements should contain the default value + assert(floats[0] == 6.0f); + assert(floats[1] == 6.0f); + assert(floats[2] == 6.0f); + + // Update values for a few elements + floats[0] = 1.0f; + floats[1] = 2.0f; + floats[2] = 3.0f; + integers[2] = -2; + assert(floats[0] == 1.0f); + assert(floats[1] == 2.0f); + assert(floats[2] == 3.0f); + assert(integers[2] == -2); + + // Reset an element, and all of its properties should revert to the defaults + properties.reset(2); + assert(floats[2] == 6.0f); + assert(integers[2] == 5); + + // Erase an element, and the size should be reduced + properties.erase(1); + assert(properties.size() == 2); + assert(properties.capacity() == 100); + assert(properties.active_list().size() == 2); + assert(properties.inactive_list().size() == 98); + + // A newly emplaced element should take the empty slot + assert(properties.emplace() == 1); + assert(properties.size() == 3); + // todo: should the new element have default properties? + assert(properties.emplace() == 3); + assert(properties.size() == 4); + + // Swapping a pair of elements swaps all of their properties + properties.swap(0, 3); + assert(integers[0] == 5); + assert(floats[0] == 6.0f); + assert(integers[3] == 5); + assert(floats[3] == 1.0f); + +} + +void test_emplace_group() { + + Property_container properties; + + auto& a = properties.add_property("a", 5); + CGAL_USE(a); + // Insert a group of 100 elements + properties.emplace_group(100); + assert(properties.size() == 100); + + // Eliminate a few regions + properties.erase(3); + assert(properties.is_erased(3)); + assert(properties.size() == 99); + for (int i = 20; i < 25; ++i) + properties.erase(i); + assert(properties.is_erased(23)); + assert(properties.size() == 94); + for (int i = 50; i < 80; ++i) + properties.erase(i); + assert(properties.is_erased(53)); + assert(properties.size() == 64); + + // A group of size 4 should only fit in the empty region fo size 5 + assert(properties.emplace_group(4) == 20); + assert(properties.size() == 68); + assert(properties.capacity() == 100); + + // A group of size 16 should only fit in the empty region fo size 30 + assert(properties.emplace_group(16) == 50); + assert(properties.size() == 84); + assert(properties.capacity() == 100); + + // Another group of size 16 should require the storage to expand, because the largest empty region is mostly full now + assert(properties.emplace_group(16) == 100); + assert(properties.size() == 100); + assert(properties.capacity() == 116); + +} + +void test_append() { + + // Create a pair of property containers with similar contents + Property_container properties_a, properties_b; + properties_a.add_property("ints", 1); + properties_b.add_property("ints", 2); + properties_a.add_property("floats", 3.0f); + properties_b.add_property("floats", 4.0f); + + // One container will also contain an extra property + properties_a.add_property("bools", true); + + // Add some values to both property sets + properties_a.emplace_group(10); + properties_b.emplace_group(5); + assert(properties_a.size() == 10); + assert(properties_b.size() == 5); + + // Add the second group to the end of the first + properties_a.append(properties_b); + assert(properties_a.size() == 15); + assert(properties_b.size() == 5); + + // Initialized values from the second group should appear after those of the first + assert(properties_a.get_property("ints")[5] == 1); + assert(properties_a.get_property("ints")[12] == 2); + assert(properties_a.get_property("floats")[5] == 3.0f); + assert(properties_a.get_property("floats")[12] == 4.0f); + + // Additional properties in the first group should have expanded too, and been filled with defaults + // note: the property array must be const, because non const operator[] doesn't work for vector! + assert(std::as_const(properties_a).get_property("bools")[12] == true); +} + +void test_constructors() { + + // Default constructor should have no properties + Property_container a{}; + assert(a.num_properties() == 0); + + // Copy constructor should duplicate all properties + auto& a_ints = a.add_property("ints", 0); + auto& a_floats = a.add_property("floats", 0.0f); + a.emplace_group(10); + a.get_property("ints")[3] = 1; + a.get_property("floats")[3] = 1.0f; + Property_container b{a}; + assert(b.num_properties() == a.num_properties() && b.num_properties() == 2); + assert(b.get_property("ints")[3] == a.get_property("ints")[3] && b.get_property("ints")[3] == 1); + assert(b.get_property("floats")[3] == a.get_property("floats")[3] && b.get_property("floats")[3] == 1.0f); + + // Copy-assignment operator should do effectively the same thing as the copy constructor + Property_container c; + c = a; + assert(c.num_properties() == a.num_properties() && c.num_properties() == 2); + assert(c.get_property("ints")[3] == a.get_property("ints")[3] && c.get_property("ints")[3] == 1); + assert(c.get_property("floats")[3] == a.get_property("floats")[3] && c.get_property("floats")[3] == 1.0f); + + // Copied property containers should not be synced with the original + a.add_property("more_ints", 2); + assert(a.num_properties() == 3); + assert(b.num_properties() == 2); + assert(c.num_properties() == 2); + a.get_property("ints")[4] = 2; + assert(a.get_property("ints")[4] == 2); + assert(b.get_property("ints")[4] == 0); + assert(c.get_property("ints")[4] == 0); + + // Copy assignment should not invalidate previously obtained array references, + // but it should update their values + auto &b_ints = b.get_property("ints"); + auto &b_floats = b.get_property("floats"); + assert(b_ints[4] == 0); + b = a; + assert(b.num_properties() == 3); + assert(b_ints[4] == 2); + + // Move assignment shouldn't invalidate references either + Property_container d{c}; + auto &d_ints = d.get_property("ints"); + assert(d_ints[4] == 0); + d = std::move(a); + assert(d.num_properties() == 3); + assert(d_ints[4] == 2); + + // Moved-from should be empty + // All properties are preserved, though + assert(a.num_properties() == 3); + assert(a.size() == 0); + assert(a_ints.capacity() == 0); + assert(a_floats.capacity() == 0); + + // Move constructor should behave like move assignment + Property_container e{std::move(b)}; + assert(e.num_properties() == 3); + assert(b.num_properties() == 3); + assert(b.size() == 0); + assert(b_ints.capacity() == 0); + assert(b_floats.capacity() == 0); +} + + +int main() { + + test_property_creation(); + test_element_access(); + test_emplace_group(); + test_append(); + test_constructors(); + + return 0; +} diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h index 88089e89e01a..6f0ab44c7ac8 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC.h @@ -253,12 +253,12 @@ class Efficient_RANSAC { /*! Retrieves the point property map. */ - const Point_map &point_map() const { return m_point_pmap; } + const Point_map &point_map() const { return *m_point_pmap; } /*! Retrieves the normal property map. */ - const Normal_map &normal() const { return m_normal_pmap; } + const Normal_map &normal() const { return *m_normal_pmap; } Input_iterator input_iterator_first() const { return m_input_iterator_first; @@ -361,13 +361,13 @@ class Efficient_RANSAC { m_direct_octrees[s] = new Direct_octree( m_traits, last + 1, last + subsetSize + 1, - m_point_pmap, + *m_point_pmap, remainingPoints - subsetSize); } else m_direct_octrees[0] = new Direct_octree( m_traits, m_input_iterator_first, m_input_iterator_first + (subsetSize), - m_point_pmap, + *m_point_pmap, 0); m_available_octree_sizes[s] = subsetSize; @@ -378,7 +378,7 @@ class Efficient_RANSAC { m_global_octree = new Indexed_octree( m_traits, m_input_iterator_first, m_input_iterator_beyond, - m_point_pmap + *m_point_pmap ); m_global_octree->refine(m_options.cluster_epsilon); @@ -496,7 +496,7 @@ class Efficient_RANSAC { } // Use bounding box diagonal as reference for default values - Bbox_3 bbox = m_global_octree->boundingBox(); + auto bbox = m_global_octree->boundingBox(); FT bbox_diagonal = (FT) CGAL::sqrt( (bbox.xmax() - bbox.xmin()) * (bbox.xmax() - bbox.xmin()) + (bbox.ymax() - bbox.ymin()) * (bbox.ymax() - bbox.ymin()) @@ -574,14 +574,14 @@ class Efficient_RANSAC { static_cast(m_num_available_points)); while (m_shape_index[first_sample] != -1); - done = drawSamplesFromCellContainingPoint - (m_global_octree, - get(m_point_pmap, - *(m_input_iterator_first + first_sample)), - select_random_octree_level(), - indices, - m_shape_index, - m_required_samples); + done = drawSamplesFromCellContainingPoint( + m_global_octree, + get(*m_point_pmap, *(m_input_iterator_first + first_sample)), + select_random_octree_level(), + indices, + m_shape_index, + m_required_samples + ); if (callback && !callback(num_invalid / double(m_num_total_points))) { clear(num_invalid, candidates); @@ -605,8 +605,8 @@ class Efficient_RANSAC { p->compute(indices, m_input_iterator_first, m_traits, - m_point_pmap, - m_normal_pmap, + *m_point_pmap, + *m_normal_pmap, m_options.epsilon, m_options.normal_threshold); @@ -1086,7 +1086,7 @@ class Efficient_RANSAC { Cell cell = stack.top(); stack.pop(); - FT width = octree->width() / (1 << (cell.depth())); + FT width = octree->width() / (1 << (octree->depth(cell))); FT diag = CGAL::sqrt(FT(3) * width * width) + epsilon; @@ -1097,10 +1097,10 @@ class Efficient_RANSAC { // differ between full or partial overlap? // if full overlap further traversal of this branch is not necessary - if (cell.is_leaf()) { + if (octree->is_leaf(cell)) { std::vector indices; - indices.reserve(cell.size()); - for (std::size_t i = 0; i < cell.size(); i++) { + indices.reserve(octree->points(cell).size()); + for (std::size_t i = 0; i < octree->points(cell).size(); i++) { if (shapeIndex[octree->index(cell, i)] == -1) { indices.push_back(octree->index(cell, i)); } @@ -1111,10 +1111,10 @@ class Efficient_RANSAC { indices); } else { - if (!cell.is_leaf()) { + if (!octree->is_leaf(cell)) { for (std::size_t i = 0; i < 8; i++) { - if (!cell[i].empty()) - stack.push(cell[i]); + if (octree->points(octree->child(cell, i)).size() != 0) + stack.push(octree->child(cell, i)); } } } @@ -1129,30 +1129,11 @@ class Efficient_RANSAC { const typename Octree::Node node_containing_point(const Octree *octree, const Point &p, std::size_t level) { // Find the node containing the point - typename Octree::Node cur = octree->root(); - while (!cur.is_null() && cur.depth() < level) { + typename Octree::Node n = octree->locate(p); + while (octree->depth(n) > level) + n = octree->parent(n); - // If cur is a leaf node, its child is null - if (cur.is_leaf()) - return typename Octree::Node(); - - // If that child is empty, return null - if (cur.empty()) - return typename Octree::Node(); - - // Determine the coordinate of the child - Point center = octree->barycenter(cur); - std::bitset<3> coordinate; - coordinate[0] = center.x() <= p.x(); - coordinate[1] = center.y() <= p.y(); - coordinate[2] = center.z() <= p.z(); - - // Otherwise, return the correct child of cur - cur = cur[coordinate.to_ulong()]; - - } - - return cur; + return n; } template @@ -1167,13 +1148,9 @@ class Efficient_RANSAC { const Cell cur = node_containing_point(octree, p, level); - // Stop if the node we need doesn't exist - if (cur.is_null()) - return false; - // Count point indices that map to -1 in the shape index std::size_t enough = 0; - for (auto j : cur) { + for (const auto j : octree->points(cur)) { if (shapeIndex[j] == -1) enough++; if (enough >= requiredSamples) @@ -1186,7 +1163,7 @@ class Efficient_RANSAC { do { std::size_t p = CGAL::get_default_random(). - uniform_int(0, cur.size() - 1); + uniform_int(0, octree->points(cur).size() - 1); std::size_t j = octree->index(cur, p); if (shapeIndex[j] == -1) @@ -1225,8 +1202,8 @@ class Efficient_RANSAC { // iterators of input data bool m_valid_iterators; Input_iterator m_input_iterator_first, m_input_iterator_beyond; - Point_map m_point_pmap; - Normal_map m_normal_pmap; + std::optional m_point_pmap; + std::optional m_normal_pmap; }; diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC_traits.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC_traits.h index 37df7a5bcc04..aa14af527c15 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC_traits.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Efficient_RANSAC_traits.h @@ -39,6 +39,8 @@ namespace CGAL { class InputPointMap, class InputNormalMap> struct Efficient_RANSAC_traits { + /// + typedef Gt GeomTraits; /// typedef typename Gt::FT FT; /// diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Octree.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Octree.h index 8e3e261e9908..acb782b70a63 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Octree.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Octree.h @@ -32,7 +32,7 @@ namespace CGAL { namespace Shape_detection { // Forward declaration needed for automatic traits detection without -// including the deprecated header itself… +// including the deprecated header itself template struct Shape_detection_traits; @@ -56,8 +56,9 @@ class RANSAC_octree { typedef std::vector Input_range; typedef Random_index_access_property_map Indexed_point_map; - typedef CGAL::Octree::type, - Input_range, Indexed_point_map> Octree; + typedef Orthtree_traits_point::type, Input_range, Indexed_point_map, false, 3> OTraits; + + typedef CGAL::Orthtree Octree; Traits m_traits; Input_range m_input_range; @@ -70,7 +71,8 @@ class RANSAC_octree { public: - typedef typename Octree::Node Node; + typedef typename Octree::Node_index Node; + typedef typename OTraits::Node_data Node_data; RANSAC_octree(const Traits &traits, Input_iterator begin, @@ -81,18 +83,26 @@ class RANSAC_octree { m_input_range(boost::counting_iterator(0), boost::counting_iterator(end - begin)), m_index_map(begin, point_map), - m_octree(m_input_range, m_index_map, 1.0), + m_octree(OTraits(m_input_range, m_index_map)), m_bBox (bbox_3(make_transform_iterator_from_property_map(begin, point_map), make_transform_iterator_from_property_map(end, point_map))), m_offset(offset) {} std::size_t index (Node node, std::size_t i) const { - return m_offset + *(node.begin() + i); + return m_offset + *(m_octree.data(node).begin() + i); + } + + std::size_t depth(const Node& node) const { + return m_octree.depth(node); + } + + bool is_leaf(const Node& node) const { + return m_octree.is_leaf(node); } std::size_t size() const { - return m_octree.root().size(); + return m_input_range.size(); } std::size_t maxLevel() const { @@ -127,17 +137,27 @@ class RANSAC_octree { return m_width; } + Node child(const Node& node, std::size_t i) const { + return m_octree.child(node, i); + } + + Node parent(const Node& node) const { + return m_octree.parent(node); + } + Node locate(const typename Traits::Point_3 &p) const { return m_octree.locate(p); } Node root() const { return m_octree.root(); } + Node_data points(const Node& n) const { return m_octree.data(n); } + typename Traits::Point_3 barycenter(const Node &node) const { return m_octree.barycenter(node); } - Bbox_3 boundingBox() const { + typename Traits::GeomTraits::Iso_cuboid_3 boundingBox() const { return m_octree.bbox(m_octree.root()); } }; diff --git a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Shape_base.h b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Shape_base.h index 8a23794aa420..99fef90cc240 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Shape_base.h +++ b/Shape_detection/include/CGAL/Shape_detection/Efficient_RANSAC/Shape_base.h @@ -374,7 +374,7 @@ namespace CGAL { */ typename boost::property_traits< typename Traits::Point_map >::reference point(std::size_t i) const { - return get(this->m_point_pmap, *(this->m_first + i)); + return get(this->m_point_pmap.value(), *(this->m_first + i)); } /*! @@ -382,7 +382,7 @@ namespace CGAL { */ typename boost::property_traits< typename Traits::Normal_map >::reference normal(std::size_t i) const { - return get(this->m_normal_pmap, *(this->m_first + i)); + return get(this->m_normal_pmap.value(), *(this->m_first + i)); } /*! @@ -694,8 +694,8 @@ namespace CGAL { Input_iterator m_first; Traits m_traits; - Point_map m_point_pmap; - Normal_map m_normal_pmap; + std::optional m_point_pmap; + std::optional m_normal_pmap; /// \endcond }; } diff --git a/Shape_detection/package_info/Shape_detection/dependencies b/Shape_detection/package_info/Shape_detection/dependencies index 63e6337e4fbf..7758a93501ea 100644 --- a/Shape_detection/package_info/Shape_detection/dependencies +++ b/Shape_detection/package_info/Shape_detection/dependencies @@ -6,6 +6,7 @@ Circulator Distance_2 Distance_3 Filtered_kernel +Hash_map Homogeneous_kernel Installation Intersections_2 @@ -16,6 +17,7 @@ Kernel_d Modular_arithmetic Number_types Orthtree +Point_set_2 Principal_component_analysis Principal_component_analysis_LGPL Profiling_tools @@ -25,4 +27,7 @@ STL_Extension Shape_detection Solver_interface Spatial_searching +Spatial_sorting Stream_support +TDS_2 +Triangulation_2 \ No newline at end of file diff --git a/Shape_regularization/package_info/Shape_regularization/dependencies b/Shape_regularization/package_info/Shape_regularization/dependencies index 5cf917f32ed4..b07b869d7157 100644 --- a/Shape_regularization/package_info/Shape_regularization/dependencies +++ b/Shape_regularization/package_info/Shape_regularization/dependencies @@ -1,5 +1,6 @@ Algebraic_foundations BGL +Cartesian_kernel Circulator Distance_2 Distance_3 @@ -14,6 +15,7 @@ Kernel_d Modular_arithmetic Number_types Orthtree +Point_set_2 Principal_component_analysis_LGPL Profiling_tools Property_map