EBGeometry_CSGOps.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_CSGOPS_HPP
#define EBGEOMETRY_CSGOPS_HPP
 
// Std includes
#include <cmath>
#include <cstdint>
#include <type_traits>
 
// Our includes
#include "EBGeometry_Affine.hpp"
#include "EBGeometry_Blend.hpp"
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_GPU.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
namespace CSGOps {
 
inline constexpr uint32_t NoTag = 0xFFFFFFFFU;
 
template <class T>
struct Tagged
{
  T m_value;
 
  uint32_t m_tag;
 
  Vec3T<T> m_point;
};
 
template <class V>
inline constexpr bool isTagged = false;
 
template <class T>
inline constexpr bool isTagged<Tagged<T>> = true;
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline T
valueOf(const T a_value) noexcept
{
  return a_value;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline T
valueOf(const Tagged<T>& a_value) noexcept
{
  return a_value.m_value;
}
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline T
withValue(const T /*a_like*/, const T a_value) noexcept
{
  return a_value;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline Tagged<T>
withValue(const Tagged<T>& a_like, const T a_value) noexcept
{
  return Tagged<T>{a_value, a_like.m_tag, a_like.m_point};
}
 
template <class V, class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
leafValue(const T a_value, [[maybe_unused]] const uint32_t a_tag, [[maybe_unused]] const Vec3T<T>& a_point) noexcept
{
  if constexpr (isTagged<V>) {
    return V{a_value, a_tag, a_point};
  }
  else {
    return a_value;
  }
}
 
template <class V>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
negated(const V& a_value) noexcept
{
  if constexpr (isTagged<V>) {
    return withValue(a_value, -a_value.m_value);
  }
  else {
    return -a_value;
  }
}
 
template <class V>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
magnitude(const V& a_value) noexcept
{
  if constexpr (isTagged<V>) {
    return withValue(a_value, std::abs(a_value.m_value));
  }
  else {
    return std::abs(a_value);
  }
}
 
template <class V, class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
plus(const V& a_value, const T a_constant) noexcept
{
  if constexpr (isTagged<V>) {
    return withValue(a_value, a_value.m_value + a_constant);
  }
  else {
    return a_value + a_constant;
  }
}
 
template <class V, class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
times(const V& a_value, const T a_constant) noexcept
{
  if constexpr (isTagged<V>) {
    return withValue(a_value, a_value.m_value * a_constant);
  }
  else {
    return a_value * a_constant;
  }
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline Tagged<T>
select(const bool a_second, const Tagged<T>& a_first, const Tagged<T>& a_other) noexcept
{
  return Tagged<T>{a_second ? a_other.m_value : a_first.m_value,
                   a_second ? a_other.m_tag : a_first.m_tag,
                   Vec3T<T>(a_second ? a_other.m_point[0] : a_first.m_point[0],
                            a_second ? a_other.m_point[1] : a_first.m_point[1],
                            a_second ? a_other.m_point[2] : a_first.m_point[2])};
}
 
template <class V>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline V
noLeaf() noexcept
{
  if constexpr (isTagged<V>) {
    using T = decltype(V::m_value);
 
    return V{T(0), NoTag, Vec3T<T>::zeros()};
  }
  else {
    return V(0);
  }
}
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline T
smaller(const T a_first, const T a_second) noexcept
{
  return Math::min(a_first, a_second);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline Tagged<T>
smaller(const Tagged<T>& a_first, const Tagged<T>& a_second) noexcept
{
  return select(a_second.m_value < a_first.m_value, a_first, a_second);
}
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline T
larger(const T a_first, const T a_second) noexcept
{
  return Math::max(a_first, a_second);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_HOST_ALWAYS_INLINE
inline Tagged<T>
larger(const Tagged<T>& a_first, const Tagged<T>& a_second) noexcept
{
  return select(a_first.m_value < a_second.m_value, a_first, a_second);
}
 
template <class T>
struct AffineCoefficients
{
  T m_values[12];
};
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline AffineCoefficients<T>
affineCoefficients(const Affine3T<T>& a_map) noexcept
{
  return AffineCoefficients<T>{{a_map.m_row0[0],
                                a_map.m_row0[1],
                                a_map.m_row0[2],
                                a_map.m_row1[0],
                                a_map.m_row1[1],
                                a_map.m_row1[2],
                                a_map.m_row2[0],
                                a_map.m_row2[1],
                                a_map.m_row2[2],
                                a_map.m_shift[0],
                                a_map.m_shift[1],
                                a_map.m_shift[2]}};
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
affinePoint(const T* a_coefficients, const Vec3T<T>& a_point) noexcept
{
  const T* m = a_coefficients;
 
  return Vec3T<T>(m[0] * a_point[0] + m[1] * a_point[1] + m[2] * a_point[2] + m[9],
                  m[3] * a_point[0] + m[4] * a_point[1] + m[5] * a_point[2] + m[10],
                  m[6] * a_point[0] + m[7] * a_point[1] + m[8] * a_point[2] + m[11]);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
affinePoint(const Affine3T<T>& a_map, const Vec3T<T>& a_point) noexcept
{
  const AffineCoefficients<T> coefficients = affineCoefficients(a_map);
 
  return affinePoint(coefficients.m_values, a_point);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
elongatePoint(const Vec3T<T>& a_halfLengths, const Vec3T<T>& a_point) noexcept
{
  return a_point - clamp(a_point,
                         Vec3T<T>(-a_halfLengths[0], -a_halfLengths[1], -a_halfLengths[2]),
                         Vec3T<T>(a_halfLengths[0], a_halfLengths[1], a_halfLengths[2]));
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
offset(const T a_value, const T a_offset) noexcept
{
  return a_value + (-a_offset);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
annular(const T a_value, const T a_delta) noexcept
{
  return std::abs(a_value) + (-a_delta);
}
 
template <class V>
struct Pair
{
  V m_first;
 
  V m_second;
};
 
template <class V>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Pair<V>
pairInit(const bool a_largest) noexcept
{
  if constexpr (isTagged<V>) {
    using T = decltype(V::m_value);
 
    const T empty = a_largest ? -Math::Limits<T>::infinity() : Math::Limits<T>::infinity();
    const V tagged{empty, NoTag, Vec3T<T>::zeros()};
 
    return Pair<V>{tagged, tagged};
  }
  else {
    const V empty = a_largest ? -Math::Limits<V>::infinity() : Math::Limits<V>::infinity();
 
    return Pair<V>{empty, empty};
  }
}
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
EBGEOMETRY_HOST_DEVICE
inline void
pairInsert(Pair<T>& a_pair, const T a_value, const bool a_largest) noexcept
{
  const bool beyondFirst  = a_largest ? (a_value > a_pair.m_first) : (a_value < a_pair.m_first);
  const bool beyondSecond = a_largest ? (a_value > a_pair.m_second) : (a_value < a_pair.m_second);
 
  if (beyondFirst) {
    a_pair.m_second = a_pair.m_first;
    a_pair.m_first  = a_value;
  }
  else if (beyondSecond) {
    a_pair.m_second = a_value;
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
pairInsert(Pair<Tagged<T>>& a_pair, const Tagged<T>& a_value, const bool a_largest) noexcept
{
  const T value  = a_value.m_value;
  const T first  = a_pair.m_first.m_value;
  const T second = a_pair.m_second.m_value;
 
  const bool beyondFirst  = a_largest ? (value > first) : (value < first);
  const bool beyondSecond = a_largest ? (value > second) : (value < second);
 
  if (beyondFirst) {
    a_pair.m_second = a_pair.m_first;
    a_pair.m_first  = a_value;
  }
  else if (beyondSecond) {
    a_pair.m_second = a_value;
  }
}
 
template <class T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
pairBlend(const Pair<T>& a_pair, const T a_smooth, const bool a_largest, const bool a_exponential) noexcept
{
  if (a_exponential) {
    return a_largest ? ExpMaxOp<T>{}(a_pair.m_first, a_pair.m_second, a_smooth)
                     : ExpMinOp<T>{}(a_pair.m_first, a_pair.m_second, a_smooth);
  }
 
  return a_largest ? SmoothMaxOp<T>{}(a_pair.m_first, a_pair.m_second, a_smooth)
                   : SmoothMinOp<T>{}(a_pair.m_first, a_pair.m_second, a_smooth);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Tagged<T>
pairBlend(const Pair<Tagged<T>>& a_pair, const T a_smooth, const bool a_largest, const bool a_exponential) noexcept
{
  const Pair<T> values{a_pair.m_first.m_value, a_pair.m_second.m_value};
 
  return withValue(a_pair.m_first, pairBlend(values, a_smooth, a_largest, a_exponential));
}
 
template <class T>
struct BlurWeights
{
  T m_values[5];
};
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST
inline BlurWeights<T>
blurWeights(const T a_distance, const T a_alpha)
{
  const T beta = (T(1) - a_alpha) / T(2);
 
  return BlurWeights<T>{{a_distance,
                         std::pow(a_alpha, T(3)) * std::pow(beta, T(0)),
                         std::pow(a_alpha, T(2)) * std::pow(beta, T(1)),
                         std::pow(a_alpha, T(1)) * std::pow(beta, T(2)),
                         std::pow(a_alpha, T(0)) * std::pow(beta, T(3))}};
}
 
template <class T, class F>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline auto
blur(const T* a_weights, const Vec3T<T>& a_point, const F& a_operand) noexcept
{
  using V = decltype(a_operand(a_point));
 
  const T  d = a_weights[0];
  const T* w = a_weights;
 
  // The 27 taps p + d (i, j, k), in this order; a tap m = |i| + |j| + |k| steps from the center has weight
  // w[1 + m]. One loop that computes each tap's point itself: with the point computed in a callable
  // called once per tap, nvcc 12.5 placed some taps at the wrong points.
  T value = T(0);
 
  if constexpr (isTagged<V>) {
    // The same sum, carrying the center tap's tag (the heaviest tap only for alpha >= 1/3).
    V center = noLeaf<V>();
 
    for (int t = 0; t < 27; t++) {
      const int i = t / 9 - 1;
      const int j = (t / 3) % 3 - 1;
      const int k = t % 3 - 1;
      const int m = i * i + j * j + k * k;
 
      const V tapValue = a_operand(a_point + Vec3T<T>(T(i) * d, T(j) * d, T(k) * d));
 
      value += w[1 + m] * tapValue.m_value;
 
      if (m == 0) {
        center = tapValue;
      }
    }
 
    return withValue(center, value);
  }
  else {
    // Kept exactly as written for the plain value: holding each tap's value in a variable first made the
    // device compiler give a kernel that evaluates a tape in a loop 255 registers instead of 168.
    for (int t = 0; t < 27; t++) {
      const int i = t / 9 - 1;
      const int j = (t / 3) % 3 - 1;
      const int k = t % 3 - 1;
      const int m = i * i + j * j + k * k;
 
      const Vec3T<T> tap = a_point + Vec3T<T>(T(i) * d, T(j) * d, T(k) * d);
 
      value += w[1 + m] * a_operand(tap);
    }
 
    return value;
  }
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
mollifyOffset(const T a_radius, const T a_spacing, const uint32_t a_i, const uint32_t a_j, const uint32_t a_k) noexcept
{
  return Vec3T<T>(-a_radius + T(a_i) * a_spacing, -a_radius + T(a_j) * a_spacing, -a_radius + T(a_k) * a_spacing);
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
mollifyWeight(const Vec3T<T>& a_offset, const T a_radius) noexcept
{
  const T t = Math::max(T(0), T(1) - a_offset.length2() / (T(4) * a_radius * a_radius));
 
  return t * t;
}
 
template <class T>
struct MollifyGrid
{
  T m_radius;
 
  T m_spacing;
 
  T m_sum;
 
  uint32_t m_samples;
};
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST
inline MollifyGrid<T>
mollifyGrid(const T a_radius, const uint32_t a_samples)
{
  MollifyGrid<T> grid{a_radius, T(2) * a_radius / T(a_samples - 1), T(0), a_samples};
 
  for (uint32_t i = 0; i < a_samples; i++) {
    for (uint32_t j = 0; j < a_samples; j++) {
      for (uint32_t k = 0; k < a_samples; k++) {
        grid.m_sum += mollifyWeight(mollifyOffset(a_radius, grid.m_spacing, i, j, k), a_radius);
      }
    }
  }
 
  return grid;
}
 
template <class T, class F>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline auto
mollify(const MollifyGrid<T>& a_grid, const Vec3T<T>& a_point, const F& a_operand) noexcept
{
  using V = decltype(a_operand(a_point));
 
  T value = T(0);
 
  if constexpr (isTagged<V>) {
    // The same sum, carrying the tag of the tap with the largest weight (the first, if several have it).
    V heaviest       = noLeaf<V>();
    T heaviestWeight = -T(1);
 
    for (uint32_t i = 0; i < a_grid.m_samples; i++) {
      for (uint32_t j = 0; j < a_grid.m_samples; j++) {
        for (uint32_t k = 0; k < a_grid.m_samples; k++) {
          const Vec3T<T> o = mollifyOffset(a_grid.m_radius, a_grid.m_spacing, i, j, k);
          const V        f = a_operand(a_point - o);
          const T        w = mollifyWeight(o, a_grid.m_radius) / a_grid.m_sum;
 
          value += f.m_value * w;
 
          if (w > heaviestWeight) {
            heaviest       = f;
            heaviestWeight = w;
          }
        }
      }
    }
 
    return withValue(heaviest, value);
  }
  else {
    // Kept exactly as written for the plain value (see blur()).
    for (uint32_t i = 0; i < a_grid.m_samples; i++) {
      for (uint32_t j = 0; j < a_grid.m_samples; j++) {
        for (uint32_t k = 0; k < a_grid.m_samples; k++) {
          const Vec3T<T> o = mollifyOffset(a_grid.m_radius, a_grid.m_spacing, i, j, k);
 
          value += a_operand(a_point - o) * (mollifyWeight(o, a_grid.m_radius) / a_grid.m_sum);
        }
      }
    }
 
    return value;
  }
}
 
template <class T>
struct Repetition
{
  Vec3T<T> m_period;
 
  Vec3T<T> m_low;
 
  Vec3T<T> m_high;
 
  Vec3T<T> m_boxLow;
 
  Vec3T<T> m_boxHigh;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  period(const size_t a_dir) const noexcept
  {
    return m_period[a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  low(const size_t a_dir) const noexcept
  {
    return m_low[a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  high(const size_t a_dir) const noexcept
  {
    return m_high[a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  boxLow(const size_t a_dir) const noexcept
  {
    return m_boxLow[a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  boxHigh(const size_t a_dir) const noexcept
  {
    return m_boxHigh[a_dir];
  }
};
 
template <class T>
struct RepetitionView
{
  const T* m_params;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  period(const size_t a_dir) const noexcept
  {
    return m_params[a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  low(const size_t a_dir) const noexcept
  {
    return m_params[3 + a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  high(const size_t a_dir) const noexcept
  {
    return m_params[6 + a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  boxLow(const size_t a_dir) const noexcept
  {
    return m_params[9 + a_dir];
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  boxHigh(const size_t a_dir) const noexcept
  {
    return m_params[12 + a_dir];
  }
};
 
template <class T, class R, class F>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline auto
repetition(const R& a_repetition, const Vec3T<T>& a_point, const F& a_object) noexcept
{
  using V = decltype(a_object(a_point));
 
  // The tile the point folds into.
  int home[3];
 
  for (size_t dir = 0; dir < 3; dir++) {
    home[dir] = static_cast<int>(std::round(
      Math::min(Math::max(a_point[dir] / a_repetition.period(dir), -a_repetition.low(dir)), a_repetition.high(dir))));
  }
 
  V best = a_object(Vec3T<T>(a_point[0] - a_repetition.period(0) * T(home[0]),
                             a_point[1] - a_repetition.period(1) * T(home[1]),
                             a_point[2] - a_repetition.period(2) * T(home[2])));
 
  // The tiles whose copy's box comes within max(0, best) of the point, along each axis.
  int first[3];
  int last[3];
 
  for (size_t dir = 0; dir < 3; dir++) {
    const T reach = Math::max(valueOf(best), T(0));
    const T lo    = std::ceil((a_point[dir] - reach - a_repetition.boxHigh(dir)) / a_repetition.period(dir));
    const T hi    = std::floor((a_point[dir] + reach - a_repetition.boxLow(dir)) / a_repetition.period(dir));
 
    first[dir] = static_cast<int>(Math::max(lo, -a_repetition.low(dir)));
    last[dir]  = static_cast<int>(Math::min(hi, a_repetition.high(dir)));
  }
 
  // Nested loops over i, j and k kept their counters alive across the calls into the object, which in a
  // Tape's sub-program raised every function of a device kernel to 255 registers (from 168); one counter,
  // in the same order (k fastest), keeps the kernel's.
  const int64_t ni    = last[0] - first[0] + 1;
  const int64_t nj    = last[1] - first[1] + 1;
  const int64_t nk    = last[2] - first[2] + 1;
  const int64_t count = (ni > 0 && nj > 0 && nk > 0) ? ni * nj * nk : 0;
 
  for (int64_t n = 0; n < count; n++) {
    const int i = first[0] + static_cast<int>(n / (nj * nk));
    const int j = first[1] + static_cast<int>((n / nk) % nj);
    const int k = first[2] + static_cast<int>(n % nk);
 
    if (i == home[0] && j == home[1] && k == home[2]) {
      continue;
    }
 
    const int tile[3] = {i, j, k};
    const T   bound   = Math::max(valueOf(best), T(0));
 
    T dist2 = T(0);
 
    for (size_t dir = 0; dir < 3; dir++) {
      const T shift = a_repetition.period(dir) * T(tile[dir]);
      const T below = (a_repetition.boxLow(dir) + shift) - a_point[dir];
      const T above = a_point[dir] - (a_repetition.boxHigh(dir) + shift);
      const T d     = Math::max(Math::max(below, above), T(0));
 
      dist2 += d * d;
    }
 
    if (dist2 <= bound * bound) {
      best = smaller(best,
                     a_object(Vec3T<T>(a_point[0] - a_repetition.period(0) * T(i),
                                       a_point[1] - a_repetition.period(1) * T(j),
                                       a_point[2] - a_repetition.period(2) * T(k))));
    }
  }
 
  return best;
}
 
} // namespace CSGOps
 
} // namespace EBGeometry
 
#endif