💿🐜 Antkeeper source code https://antkeeper.com
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

577 lines
16 KiB

/*
* Copyright (C) 2021 Christopher J. Howard
*
* This file is part of Antkeeper source code.
*
* Antkeeper source code is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Antkeeper source code is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Antkeeper source code. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef ANTKEEPER_GEOM_HYPEROCTREE_HPP
#define ANTKEEPER_GEOM_HYPEROCTREE_HPP
#include "math/compile.hpp"
#include <algorithm>
#include <array>
#include <bit>
#include <concepts>
#include <cstddef>
#include <set>
#include <type_traits>
#include <unordered_set>
namespace geom {
/// Orders in which hyperoctree nodes can be stored and iterated.
enum class hyperoctree_order
{
/// Hyperoctree nodes are unordered, potentially resulting in faster insertions through the internal use of `std::unordered_set` rather than `std::set`.
unordered,
/// Hyperoctree nodes are stored and iterated in depth-first preorder.
dfs_pre,
/// Hyperoctree nodes are stored and iterated in breadth-first order.
bfs
};
/// @private
template <std::unsigned_integral T>
using hyperoctree_dfs_pre_compare = std::less<T>;
/// @private
template <std::unsigned_integral T, std::size_t DepthBits>
struct hyperoctree_bfs_compare
{
constexpr bool operator()(const T& lhs, const T& rhs) const noexcept
{
return std::rotr(lhs, DepthBits) < std::rotr(rhs, DepthBits);
}
};
/// @private
template <hyperoctree_order Order, std::unsigned_integral T, std::size_t DepthBits>
struct hyperoctree_container {};
/// @private
template <std::unsigned_integral T, std::size_t DepthBits>
struct hyperoctree_container<hyperoctree_order::unordered, T, DepthBits>
{
typedef std::unordered_set<T> type;
};
/// @private
template <std::unsigned_integral T, std::size_t DepthBits>
struct hyperoctree_container<hyperoctree_order::dfs_pre, T, DepthBits>
{
typedef std::set<T, hyperoctree_dfs_pre_compare<T>> type;
};
/// @private
template <std::unsigned_integral T, std::size_t DepthBits>
struct hyperoctree_container<hyperoctree_order::bfs, T, DepthBits>
{
typedef std::set<T, hyperoctree_bfs_compare<T, DepthBits>> type;
};
/**
* Hashed linear hyperoctree.
*
* @tparam T Unsigned integral node identifier type.
* @tparam N Number of dimensions.
* @tparam Order Order in which nodes are stored and iterated.
*
* @see http://codervil.blogspot.com/2015/10/octree-node-identifiers.html
* @see https://geidav.wordpress.com/2014/08/18/advanced-octrees-2-node-representations/
*/
template <std::unsigned_integral T, std::size_t N, hyperoctree_order Order = hyperoctree_order::dfs_pre>
class hyperoctree
{
private:
/**
* Finds the maximum hyperoctree depth level from the size of the node type `T` and number of dimensions `N`.
*
* @return Maximum hyperoctree depth level.
*
* @note There is likely a more elegant formula for this. Information about the 2D and 3D cases is given below:
*
* 2D:
* 8 bit ( 1 byte) = max depth 1 ( 4 loc bits, 1 depth bits, 1 divider bit) = 6 bits
* 16 bit ( 2 byte) = max depth 5 ( 12 loc bits, 3 depth bits, 1 divider bit) = 16 bits
* 32 bit ( 4 byte) = max depth 12 ( 26 loc bits, 4 depth bits, 1 divider bit) = 31 bits
* 64 bit ( 8 byte) = max depth 28 ( 58 loc bits, 5 depth bits, 1 divider bit) = 64 bits
* 128 bit (16 byte) = max depth 59 (120 loc bits, 6 depth bits, 1 divider bit) = 127 bits
* 256 bit (32 byte) = max depth 123 (248 loc bits, 7 depth bits, 1 divider bit) = 256 bits
*
* @see https://oeis.org/A173009
*
* 3D:
* 8 bit ( 1 byte) = max depth 1 ( 6 loc bits, 1 depth bits, 1 divider bit) = 8 bits
* 16 bit ( 2 byte) = max depth 3 ( 12 loc bits, 2 depth bits, 1 divider bit) = 15 bits
* 32 bit ( 4 byte) = max depth 8 ( 27 loc bits, 4 depth bits, 1 divider bit) = 32 bits
* 64 bit ( 8 byte) = max depth 18 ( 57 loc bits, 5 depth bits, 1 divider bit) = 63 bits
* 128 bit (16 byte) = max depth 39 (120 loc bits, 6 depth bits, 1 divider bit) = 127 bits
* 256 bit (32 byte) = max depth 81 (243 loc bits, 7 depth bits, 1 divider bit) = 251 bits
*
* @see https://oeis.org/A178420
*/
static consteval std::size_t find_max_depth() noexcept
{
std::size_t max_depth = 0;
for (std::size_t i = 1; i <= sizeof(T) * 8; ++i)
{
const std::size_t location_bits = sizeof(T) * 8 - i;
max_depth = location_bits / N - 1;
const std::size_t depth_bits = static_cast<std::size_t>(std::bit_width(max_depth));
if (depth_bits + location_bits < sizeof(T) * 8)
break;
}
return static_cast<T>(max_depth);
}
public:
/// Node identifier type.
typedef T node_type;
/// Number of dimensions.
static constexpr std::size_t dimensions = N;
/// Node storage and traversal order.
static constexpr hyperoctree_order order = Order;
/// Maximum node depth level.
static constexpr node_type max_depth = find_max_depth();
/// Number of bits in the node type.
static constexpr node_type node_bits = sizeof(node_type) * 8;
/// Number of bits required to encode the depth of a node.
static constexpr node_type depth_bits = std::bit_width(max_depth);
/// Number of bits required to encode the Morton location code of a node.
static constexpr node_type location_bits = (max_depth + 1) * N;
/// Number of bits separating the depth and Morton location code in a node identifier.
static constexpr node_type divider_bits = node_bits - (depth_bits + location_bits);
/// Number of children per node.
static constexpr node_type children_per_node = math::compile::exp2<node_type>(N);
/// Number of siblings per node.
static constexpr node_type siblings_per_node = children_per_node - 1;
/// Resolution in each dimension.
static constexpr node_type resolution = math::compile::exp2<node_type>(max_depth + 1);
/// Number of nodes in a full hyperoctree.
static constexpr std::size_t max_node_count = (math::compile::pow<std::size_t>(resolution, N) - 1) / siblings_per_node;
/// Node identifier of the persistent root node.
static constexpr node_type root = 0;
/// Node container type.
typedef typename hyperoctree_container<order, node_type, depth_bits>::type container_type;
/// Iterator type.
typedef typename container_type::iterator iterator;
/// Constant iterator type.
typedef typename container_type::const_iterator const_iterator;
/// Reverse iterator type.
typedef std::conditional<order != hyperoctree_order::unordered, std::reverse_iterator<iterator>, iterator>::type reverse_iterator;
/// Constant reverse iterator type.
typedef std::conditional<order != hyperoctree_order::unordered, std::reverse_iterator<const_iterator>, const_iterator>::type const_reverse_iterator;
/// @name Nodes
/// @{
/**
* Extracts the depth of a node from its identifier.
*
* @param node Node identifier.
*
* @return Depth of the node.
*/
static constexpr inline node_type depth(node_type node) noexcept
{
constexpr node_type mask = math::compile::exp2<node_type>(depth_bits) - 1;
return node & mask;
}
/**
* Extracts the Morton location code of a node from its identifier.
*
* @param node Node identifier.
*
* @return Morton location code of the node.
*/
static constexpr inline node_type location(node_type node) noexcept
{
return node >> ((node_bits - 1) - depth(node) * N);
}
/**
* Extracts the depth and Morton location code of a node from its identifier.
*
* @param node Node identifier.
*
* @return Array containing the depth of the node, followed by the Morton location code of the node.
*/
static constexpr std::array<node_type, 2> split(node_type node) noexcept
{
const node_type depth = hyperoctree::depth(node);
const node_type location = node >> ((node_bits - 1) - depth * N);
return {depth, location};
}
/**
* Constructs an identifier for a node at the given depth and location.
*
* @param depth Depth level.
* @param location Morton location code.
*
* @return Identifier of a node at the given depth and location.
*
* @warning If @p depth exceeds `max_depth`, the returned node identifier is not valid.
*/
static constexpr inline node_type node(node_type depth, node_type location) noexcept
{
return (location << ((node_bits - 1) - depth * N)) | depth;
}
/**
* Constructs an identifier for the ancestor of a node at a given depth.
*
* @param node Node identifier.
* @param depth Absolute depth of an ancestor node.
*
* @return Identifier of the ancestor of the node at the given depth.
*
* @warning If @p depth exceeds the depth of @p node, the returned node identifier is not valid.
*/
static constexpr inline node_type ancestor(node_type node, node_type depth) noexcept
{
const node_type mask = (~node_type(0)) << ((node_bits - 1) - depth * N);
return (node & mask) | depth;
}
/**
* Constructs an identifier for the parent of a node.
*
* @param node Node identifier.
*
* @return Identifier of the parent node.
*/
static constexpr inline node_type parent(node_type node) noexcept
{
return ancestor(node, depth(node) - 1);
}
/**
* Constructs an identifier for the nth sibling of a node.
*
* @param node Node identifier.
* @param n Offset to a sibling, automatically wrapped to `[0, siblings_per_node]`.
*
* @return Identifier of the nth sibling node.
*/
static constexpr node_type sibling(node_type node, node_type n) noexcept
{
constexpr node_type mask = (1 << N) - 1;
const auto [depth, location] = split(node);
const node_type sibling_location = (location & (~mask)) | ((location + n) & mask);
return hyperoctree::node(depth, sibling_location);
}
/**
* Constructs an identifier for the nth child of a node.
*
* @param node Node identifier.
* @param n Offset to a sibling of the first child node, automatically wrapped to `[0, siblings_per_node]`.
*
* @return Identifier of the nth child node.
*/
static constexpr inline node_type child(node_type node, T n) noexcept
{
return sibling(node + 1, n);
}
/**
* Constructs an identifier for first common ancestor of two nodes
*
* @param a Identifier of the first node.
* @param b Identifier of the second node.
*
* @return Identifier of the first common ancestor of the two nodes.
*/
static constexpr node_type common_ancestor(node_type a, node_type b) noexcept
{
const node_type bits = std::min<node_type>(depth(a), depth(b)) * N;
const node_type marker = (node_type(1) << (node_bits - 1)) >> bits;
const node_type depth = node_type(std::countl_zero((a ^ b) | marker) / N);
return ancestor(a, depth);
}
/// @}
/// Constructs a hyperoctree with a single root node.
hyperoctree():
nodes({root})
{}
/// @name Iterators
/// @{
/**
* Returns an iterator to the first node, in the traversal order specified by hyperoctree::order.
*
* @note Node identifiers cannot be modified through iterators.
*/
/// @{
inline iterator begin() noexcept
{
return nodes.begin();
}
inline const_iterator begin() const noexcept
{
return nodes.begin();
}
inline const_iterator cbegin() const noexcept
{
return nodes.cbegin();
}
/// @}
/**
* Returns an iterator to the node following the last node, in the traversal order specified by hyperoctree::order.
*
* @note Node identifiers cannot be modified through iterators.
*/
/// @{
inline iterator end() noexcept
{
return nodes.end();
}
inline const_iterator end() const noexcept
{
return nodes.end();
}
inline const_iterator cend() const noexcept
{
return nodes.cend();
}
/// @}
/**
* Returns a reverse iterator to the first node of the revered hyperoctree, in the traversal order specified by hyperoctree::order.
*
* @note Node identifiers cannot be modified through iterators.
* @note If the hyperoctree is unordered, reverse iteration and forward iteration will be identical.
*/
/// @{
inline reverse_iterator rbegin() noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.rbegin();
else
return nodes.begin();
}
inline const_reverse_iterator rbegin() const noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.rbegin();
else
return nodes.begin();
}
inline const_reverse_iterator crbegin() const noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.crbegin();
else
return nodes.cbegin();
}
/// @}
/**
* Returns a reverse iterator to the node following the last node of the reverse hyperoctree, in the traversal order specified by hyperoctree::order.
*
* @note Node identifiers cannot be modified through iterators.
* @note If the hyperoctree is unordered, reverse iteration and forward iteration will be identical.
*/
/// @{
inline reverse_iterator rend() noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.rend();
else
return nodes.end();
}
inline const_reverse_iterator rend() const noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.rend();
else
return nodes.end();
}
inline const_reverse_iterator crend() const noexcept
{
if constexpr (order != hyperoctree_order::unordered)
return nodes.crend();
else
return nodes.cend();
}
/// @}
/// @}
/// @name Capacity
/// @{
/**
* Checks if the hyperoctree has no nodes.
*
* @return `true` if the hyperoctree is empty, `false` otherwise.
*
* @note This function should always return `false`, as the root node is persistent.
*/
inline bool empty() const noexcept
{
return nodes.empty();
}
/**
* Checks if the hyperoctree is full.
*
* @return `true` if the hyperoctree is full, `false` otherwise.
*/
inline bool full() const noexcept
{
return size() == max_size();
}
/**
* Returns the number of nodes in the hyperoctree.
*
* @return Number of nodes in the hyperoctree.
*
* @note Hyperoctree size will always be greater than or equal to one, as the root node is persistent.
*/
inline std::size_t size() const noexcept
{
return nodes.size();
}
/// Returns the total number of nodes the hyperoctree is capable of containing.
static consteval std::size_t max_size() noexcept
{
return max_node_count;
}
/// @}
/// @name Modifiers
/// @{
/**
* Erases all nodes except the root node, which is persistent.
*/
inline void clear()
{
nodes = {root};
}
/**
* Inserts a node and its siblings into the hyperoctree, creating ancestors as necessary.
*
* @param node Node to insert.
*
* @note The root node is persistent and does not need to be inserted.
*/
void insert(node_type node)
{
if (contains(node))
return;
// Insert node
nodes.emplace(node);
// Insert node siblings
for (node_type i = 1; i < children_per_node; ++i)
nodes.emplace(sibling(node, i));
// Insert node ancestors
insert(parent(node));
}
/**
* Erases a node, along with its descendants, siblings, and descendants of siblings.
*
* @param node Identifier of the node to erase.
*
* @note The root node is persistent and cannot be erased.
*/
void erase(node_type node)
{
if (node == root || !contains(node))
return;
// Erase node and its descendants
nodes.erase(node);
erase(child(node, 0));
// Erase node siblings
for (node_type i = 0; i < siblings_per_node; ++i)
{
node = sibling(node, 1);
// Erase sibling and its descendants
nodes.erase(node);
erase(child(node, 0));
}
}
/// @}
/// @name Lookup
/// @{
/**
* Checks if a node is contained within the hyperoctree.
*
* @param node Identifier of the node to check for.
*
* @return `true` if the hyperoctree contains the node, `false` otherwise.
*/
inline bool contains(node_type node) const
{
return nodes.contains(node);
}
/**
* Checks if a node has no children.
*
* @param node Node identififer.
*
* @return `true` if the node has no children, and `false` otherwise.
*/
inline bool is_leaf(node_type node) const
{
return !contains(child(node, 0));
}
/// @}
private:
container_type nodes;
};
/// Hyperoctree with unordered node storage and traversal.
template <std::unsigned_integral T, std::size_t N>
using unordered_hyperoctree = hyperoctree<T, N, hyperoctree_order::unordered>;
} // namespace geom
#endif // ANTKEEPER_GEOM_HYPEROCTREE_HPP