Morton.h

Go to the documentation of this file.

 
#ifndef PBAT_GEOMETRY_MORTON_H
#define PBAT_GEOMETRY_MORTON_H
 
#include "pbat/HostDevice.h"
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <type_traits>
 
namespace pbat {
namespace geometry {
 
using MortonCodeType = std::uint32_t; 

PBAT_HOST_DEVICE inline MortonCodeType ExpandBits(MortonCodeType v)
{
    v = (v * 0x00010001u) & 0xFF0000FFu;
    v = (v * 0x00000101u) & 0x0F00F00Fu;
    v = (v * 0x00000011u) & 0xC30C30C3u;
    v = (v * 0x00000005u) & 0x49249249u;
    return v;
}
 
namespace detail {
 
template <class Point>
concept CMorton3dPoint = requires(Point p)
{
    {
        p[0]
    } -> std::convertible_to<float>;
    {
        p[1]
    } -> std::convertible_to<float>;
    {
        p[2]
    } -> std::convertible_to<float>;
};
 
} // namespace detail

template <detail::CMorton3dPoint Point>
[[maybe_unused]] PBAT_HOST_DEVICE inline MortonCodeType Morton3D(Point x)
{
    using namespace std;
    using ScalarType  = std::remove_cvref_t<decltype(x[0])>;
    MortonCodeType xx = ExpandBits(static_cast<MortonCodeType>(
        min(max(x[0] * ScalarType(1024), ScalarType(0)), ScalarType(1023))));
    MortonCodeType yy = ExpandBits(static_cast<MortonCodeType>(
        min(max(x[1] * ScalarType(1024), ScalarType(0)), ScalarType(1023))));
    MortonCodeType zz = ExpandBits(static_cast<MortonCodeType>(
        min(max(x[2] * ScalarType(1024), ScalarType(0)), ScalarType(1023))));
    return xx * 4 + yy * 2 + zz;
}
 
} // namespace geometry
} // namespace pbat
 
#endif // PBAT_GEOMETRY_MORTON_H