EBGeometry_PointCloudDetail.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2026 Robert Marskar <robert.marskar@sintef.no>
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef EBGEOMETRY_POINTCLOUDDETAIL_HPP
#define EBGEOMETRY_POINTCLOUDDETAIL_HPP
 
// Std includes
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <vector>
 
// Our includes
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_GPU.hpp"
#include "EBGeometry_GPUAlgorithms.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_PointCloud.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
namespace PointCloudDetail {
 
template <class T>
inline void
requireValidCloud(const char* a_who, const std::vector<Vec3T<T>>& a_positions) noexcept
{
  EBGEOMETRY_REQUIRE(a_positions.size() < std::size_t(PointCloud::InvalidIndex),
                     "%s: too many points for 32-bit indices (%zu points, limit %zu)",
                     a_who,
                     a_positions.size(),
                     std::size_t(PointCloud::InvalidIndex) - 1);
 
  for (std::size_t i = 0; i < a_positions.size(); i++) {
    const Vec3T<T>& p = a_positions[i];
 
    EBGEOMETRY_REQUIRE(std::isfinite(p[0]) && std::isfinite(p[1]) && std::isfinite(p[2]),
                       "%s: point %zu of %zu has a non-finite coordinate",
                       a_who,
                       i,
                       a_positions.size());
  }
}
 
[[nodiscard]] EBGEOMETRY_HOST
inline uint32_t
keyCount(const uint32_t a_n) noexcept
{
  uint32_t n = (a_n == 0) ? 0 : 1;
 
  while (n < a_n) {
    n *= 2;
  }
 
  return n;
}
 
template <class T, class Source>
[[nodiscard]] EBGEOMETRY_HOST
inline BoundingVolumes::AABBT<T>
boxOf(const Source&  a_source,
      const uint32_t a_numPoints,
      const bool     a_parallel,
      Vec3T<T>*      a_scratch,
      const char*    a_who,
      const char*    a_noun = "point")
{
  if (a_numPoints == 0) {
    return BoundingVolumes::AABBT<T>(Vec3T<T>::zeros(), Vec3T<T>::zeros());
  }
 
  if (a_parallel) {
    const BoundingVolumes::AABBT<T> box = GPU::boundingBoxOf<T>(a_source, a_numPoints, a_scratch);
 
    const Vec3T<T>& lo = box.getLowCorner();
    const Vec3T<T>& hi = box.getHighCorner();
 
    EBGEOMETRY_REQUIRE(std::isfinite(lo[0]) && std::isfinite(lo[1]) && std::isfinite(lo[2]) && std::isfinite(hi[0]) &&
                         std::isfinite(hi[1]) && std::isfinite(hi[2]),
                       "%s: a %s has a non-finite coordinate",
                       a_who,
                       a_noun);
 
    return box;
  }
 
  Vec3T<T> lo = +Vec3T<T>::max();
  Vec3T<T> hi = -Vec3T<T>::max();
 
  for (uint32_t i = 0; i < a_numPoints; i++) {
    const Vec3T<T> p = a_source(i);
 
    EBGEOMETRY_REQUIRE(std::isfinite(p[0]) && std::isfinite(p[1]) && std::isfinite(p[2]),
                       "%s: %s %u of %u has a non-finite coordinate",
                       a_who,
                       a_noun,
                       unsigned(i),
                       unsigned(a_numPoints));
 
    lo = min(lo, p);
    hi = max(hi, p);
  }
 
  return BoundingVolumes::AABBT<T>(lo, hi);
}
 
} // namespace PointCloudDetail
 
} // namespace EBGeometry
 
#endif