EBGeometry_Fused.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_FUSED_HPP
#define EBGEOMETRY_FUSED_HPP
 
// Std includes
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <utility>
 
// Our includes
#include "EBGeometry_Affine.hpp"
#include "EBGeometry_BVHUnion.hpp"
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_CSGOps.hpp"
#include "EBGeometry_DistanceQuality.hpp"
#include "EBGeometry_Expr.hpp"
#include "EBGeometry_GPU.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_Pool.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
namespace Fused {
 
namespace Detail {
 
template <class... Ss>
struct Pack;
 
template <>
struct Pack<>
{
};
 
template <class S, class... Rest>
struct Pack<S, Rest...>
{
  S m_head;
 
  Pack<Rest...> m_tail;
};
 
template <size_t I, class S, class... Rest>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline const auto&
get(const Pack<S, Rest...>& a_pack) noexcept
{
  if constexpr (I == 0) {
    return a_pack.m_head;
  }
  else {
    return get<I - 1>(a_pack.m_tail);
  }
}
 
template <size_t I, class S, class... Rest>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline auto&
get(Pack<S, Rest...>& a_pack) noexcept
{
  if constexpr (I == 0) {
    return a_pack.m_head;
  }
  else {
    return get<I - 1>(a_pack.m_tail);
  }
}
 
template <size_t I, class... Cs>
struct TypeAt;
 
template <class C, class... Rest>
struct TypeAt<0, C, Rest...>
{
  using type = C;
};
 
template <size_t I, class C, class... Rest>
struct TypeAt<I, C, Rest...>
{
  using type = typename TypeAt<I - 1, Rest...>::type;
};
 
template <class... Cs>
inline constexpr bool anyPoolHandle = (isPoolHandle<Cs> || ...);
 
template <class... Cs>
[[nodiscard]] constexpr DistanceQuality
worstAndBound() noexcept
{
  const DistanceQuality qualities[] = {distanceQualityOf<Cs>...};
 
  DistanceQuality q = DistanceQuality::Bound;
 
  for (const DistanceQuality c : qualities) {
    q = (static_cast<int>(c) > static_cast<int>(q)) ? c : q;
  }
 
  return q;
}
 
template <bool Handle>
struct Base
{
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline PoolLocation
  location() const noexcept
  {
    return PoolLocation{};
  }
};
 
template <>
struct Base<true> : PoolHandle
{
  Base() noexcept = default;
 
  EBGEOMETRY_HOST_DEVICE
  explicit Base(const PoolLocation& a_location) noexcept : PoolHandle(a_location)
  {}
};
 
template <class... Cs>
struct Operands
{
  const Pack<CSGDetail::Stored<Cs>...>& m_pack;
 
  PoolLocation m_location;
 
  static constexpr size_t count = sizeof...(Cs);
 
  template <size_t I, class T>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  value(const Vec3T<T>& a_point) const noexcept
  {
    using C = typename TypeAt<I, Cs...>::type;
 
    if constexpr (isPoolHandle<C>) {
      return C(get<I>(m_pack), m_location).signedDistance(a_point);
    }
    else {
      return get<I>(m_pack).signedDistance(a_point);
    }
  }
 
  template <size_t I>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline auto
  box() const
  {
    using C = typename TypeAt<I, Cs...>::type;
 
    return CSGDetail::fromStored<C>(get<I>(m_pack), m_location).computeBoundingVolume();
  }
};
 
} // namespace Detail
 
template <class Op, class... Cs>
class Node : public Detail::Base<Detail::anyPoolHandle<Cs...>>
{
public:
  static_assert(sizeof...(Cs) > 0, "Fused::Node needs an operand");
  static_assert((isShape<Cs> && ...), "Fused::Node: every operand must be a shape (see isShape)");
 
  using T = ShapePrecision<typename Detail::TypeAt<0, Cs...>::type>;
 
  static_assert((std::is_same_v<ShapePrecision<Cs>, T> && ...), "Fused::Node: the operands' precisions differ");
 
  static constexpr bool isHandle = Detail::anyPoolHandle<Cs...>;
 
  static constexpr DistanceQuality distanceQuality = Op::template quality<Cs...>();
 
  struct Layout
  {
    Op m_op;
 
    Detail::Pack<CSGDetail::Stored<Cs>...> m_operands;
  };
 
  Node() = delete;
 
  template <bool H = isHandle, std::enable_if_t<!H, int> = 0>
  EBGEOMETRY_HOST_DEVICE
  explicit Node(const Layout& a_layout) noexcept : m_layout(a_layout)
  {}
 
  template <bool H = isHandle, std::enable_if_t<H, int> = 0>
  EBGEOMETRY_HOST_DEVICE
  Node(const Layout& a_layout, const PoolLocation& a_location) noexcept
    : Detail::Base<true>(a_location), m_layout(a_layout)
  {}
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline const Layout&
  layout() const noexcept
  {
    return m_layout;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    return m_layout.m_op.evaluate(this->operands(), a_point);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  computeBoundingVolume() const
  {
    return m_layout.m_op.box(this->operands());
  }
 
  [[nodiscard]] EBGEOMETRY_HOST
  inline Node
  rebasedView(const Pool& a_pool) const noexcept;
 
  [[nodiscard]] EBGEOMETRY_HOST
  inline Node
  deepCopy(Pool& a_dstPool) const;
 
private:
  Layout m_layout;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline Detail::Operands<Cs...>
  operands() const noexcept
  {
    return Detail::Operands<Cs...>{m_layout.m_operands, this->location()};
  }
};
 
template <bool IsUnion>
struct CombineOp
{
  template <class... Cs>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<Cs...>();
  }
 
  template <class O, class T>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    return this->fold(a_operands, a_point, std::make_index_sequence<O::count - 1>{});
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline auto
  box(const O& a_operands) const
  {
    return this->boxes(a_operands, std::make_index_sequence<O::count>{});
  }
 
private:
  template <class O, class T, size_t... Is>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  fold(const O& a_operands, const Vec3T<T>& a_point, std::index_sequence<Is...>) const noexcept
  {
    T value = a_operands.template value<0>(a_point);
 
    if constexpr (IsUnion) {
      ((value = Math::min(value, a_operands.template value<Is + 1>(a_point))), ...);
    }
    else {
      ((value = Math::max(value, a_operands.template value<Is + 1>(a_point))), ...);
    }
 
    return value;
  }
 
  template <class O, size_t... Is>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline auto
  boxes(const O& a_operands, std::index_sequence<Is...>) const
  {
    using Box = decltype(a_operands.template box<0>());
    using T   = std::decay_t<decltype(a_operands.template box<0>().getLowCorner()[0])>;
 
    Box box = IsUnion ? Box() : ExprDetail::everywhere<T>();
 
    if constexpr (IsUnion) {
      ((box = box.merged(a_operands.template box<Is>())), ...);
    }
    else {
      ((box = ExprDetail::intersected(box, a_operands.template box<Is>())), ...);
    }
 
    return box;
  }
};
 
using UnionOp = CombineOp<true>;
 
using IntersectionOp = CombineOp<false>;
 
template <class T, bool IsUnion>
struct SmoothCombineOp
{
  T m_smooth;
 
  bool m_exponential;
 
  template <class... Cs>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<Cs...>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    return this->blend(a_operands, a_point, std::make_index_sequence<O::count>{});
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    const BoundingVolumes::AABBT<T> box = CombineOp<IsUnion>{}.box(a_operands);
 
    if constexpr (IsUnion) {
      return ExprDetail::grown(box, m_exponential ? m_smooth * std::log(T(2)) : T(0.25) * m_smooth);
    }
    else {
      return box;
    }
  }
 
private:
  template <class O, size_t... Is>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  blend(const O& a_operands, const Vec3T<T>& a_point, std::index_sequence<Is...>) const noexcept
  {
    CSGOps::Pair<T> pair = CSGOps::pairInit<T>(!IsUnion);
 
    (CSGOps::pairInsert(pair, a_operands.template value<Is>(a_point), !IsUnion), ...);
 
    return CSGOps::pairBlend(pair, m_smooth, !IsUnion, m_exponential);
  }
};
 
struct ComplementOp
{
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return distanceQualityOf<C>;
  }
 
  template <class O, class T>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    return -a_operands.template value<0>(a_point);
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline auto
  box(const O& a_operands) const
  {
    using T = std::decay_t<decltype(a_operands.template box<0>().getLowCorner()[0])>;
 
    return ExprDetail::everywhere<T>();
  }
};
 
template <class T>
struct AffineOp
{
  Affine3T<T> m_map;
 
  Affine3T<T> m_inverse;
 
  T m_valueScale;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return distanceQualityOf<C>;
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    const T value = a_operands.template value<0>(CSGOps::affinePoint(m_map, a_point));
 
    return (m_valueScale != T(1)) ? value * m_valueScale : value;
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    return m_inverse.transform(a_operands.template box<0>());
  }
};
 
template <class T>
struct ElongateOp
{
  Vec3T<T> m_halfLengths;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<C>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    return a_operands.template value<0>(CSGOps::elongatePoint(m_halfLengths, a_point));
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    return ExprDetail::grown(a_operands.template box<0>(), m_halfLengths);
  }
};
 
template <class T, bool IsAnnulus>
struct ShiftOp
{
  T m_amount;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return IsAnnulus ? distanceQualityOf<C> : Detail::worstAndBound<C>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    const T value = a_operands.template value<0>(a_point);
 
    if constexpr (IsAnnulus) {
      return CSGOps::annular(value, m_amount);
    }
    else {
      return CSGOps::offset(value, m_amount);
    }
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    return ExprDetail::grown(a_operands.template box<0>(), Math::max(m_amount, T(0)));
  }
};
 
template <class T>
struct RepetitionOp
{
  CSGOps::Repetition<T> m_repetition;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<C>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    const auto f = [&a_operands](const Vec3T<T>& a_x) noexcept -> T { return a_operands.template value<0>(a_x); };
 
    return CSGOps::repetition<T>(m_repetition, a_point, f);
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const;
};
 
template <class T>
struct BlurOp
{
  CSGOps::BlurWeights<T> m_weights;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<C>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    const auto f = [&a_operands](const Vec3T<T>& a_x) noexcept -> T { return a_operands.template value<0>(a_x); };
 
    return CSGOps::blur(m_weights.m_values, a_point, f);
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    return ExprDetail::grown(a_operands.template box<0>(), std::sqrt(T(3)) * m_weights.m_values[0]);
  }
};
 
template <class T>
struct MollifyOp
{
  CSGOps::MollifyGrid<T> m_grid;
 
  bool m_identity;
 
  template <class C>
  [[nodiscard]] static constexpr DistanceQuality
  quality() noexcept
  {
    return Detail::worstAndBound<C>();
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  evaluate(const O& a_operands, const Vec3T<T>& a_point) const noexcept
  {
    if (m_identity) {
      return a_operands.template value<0>(a_point);
    }
 
    const auto f = [&a_operands](const Vec3T<T>& a_x) noexcept -> T { return a_operands.template value<0>(a_x); };
 
    return CSGOps::mollify(m_grid, a_point, f);
  }
 
  template <class O>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  box(const O& a_operands) const
  {
    const BoundingVolumes::AABBT<T> box = a_operands.template box<0>();
 
    return m_identity ? box : ExprDetail::grown(box, std::sqrt(T(3)) * m_grid.m_radius);
  }
};
 
template <class... Cs>
struct OperandList
{
  Detail::Pack<Cs...> m_items;
};
 
template <class... Cs>
[[nodiscard]] EBGEOMETRY_HOST
inline OperandList<Cs...>
operands(const Cs&... a_operands);
 
template <class Op, class... Cs>
using NodeOf = Node<Op, Cs...>;
 
template <class C, class... Cs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Union(const C& a_first, const Cs&... a_others);
 
template <class C, class... Cs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Intersection(const C& a_first, const Cs&... a_others);
 
template <class A, class B, class... Bs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Difference(const A& a_first, const B& a_second, const Bs&... a_others);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<ComplementOp, C>
Complement(const C& a_operand);
 
template <class A, class B>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<SmoothCombineOp<ShapePrecision<A>, true>, A, B>
SmoothUnion(const A&                a_first,
            const B&                a_second,
            const ShapePrecision<A> a_smooth,
            const BlendKind         a_blend = BlendKind::Polynomial);
 
template <class... Cs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
SmoothUnion(const OperandList<Cs...>&                                     a_operands,
            const ShapePrecision<typename Detail::TypeAt<0, Cs...>::type> a_smooth,
            const BlendKind                                               a_blend = BlendKind::Polynomial);
 
template <class A, class B>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<SmoothCombineOp<ShapePrecision<A>, false>, A, B>
SmoothIntersection(const A&                a_first,
                   const B&                a_second,
                   const ShapePrecision<A> a_smooth,
                   const BlendKind         a_blend = BlendKind::Polynomial);
 
template <class... Cs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
SmoothIntersection(const OperandList<Cs...>&                                     a_operands,
                   const ShapePrecision<typename Detail::TypeAt<0, Cs...>::type> a_smooth,
                   const BlendKind                                               a_blend = BlendKind::Polynomial);
 
template <class A, class B>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
SmoothDifference(const A&                a_first,
                 const B&                a_second,
                 const ShapePrecision<A> a_smooth,
                 const BlendKind         a_blend = BlendKind::Polynomial);
 
template <class A, class... Bs>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
SmoothDifference(const A&                  a_first,
                 const OperandList<Bs...>& a_others,
                 const ShapePrecision<A>   a_smooth,
                 const BlendKind           a_blend = BlendKind::Polynomial);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Translate(const C& a_operand, const Vec3T<ShapePrecision<C>>& a_shift);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Rotate(const C& a_operand, const ShapePrecision<C> a_degrees, const size_t a_axis);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Rotate(const C& a_operand, const ShapePrecision<C> a_degrees, const Vec3T<ShapePrecision<C>>& a_axis);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Scale(const C& a_operand, const ShapePrecision<C> a_factor);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline auto
Reflect(const C& a_operand, const size_t a_axis);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<ShiftOp<ShapePrecision<C>, false>, C>
Offset(const C& a_operand, const ShapePrecision<C> a_offset);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<ShiftOp<ShapePrecision<C>, true>, C>
Annular(const C& a_operand, const ShapePrecision<C> a_delta);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<ElongateOp<ShapePrecision<C>>, C>
Elongate(const C& a_operand, const Vec3T<ShapePrecision<C>>& a_elongation);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<RepetitionOp<ShapePrecision<C>>, C>
FiniteRepetition(const C&                        a_operand,
                 const Vec3T<ShapePrecision<C>>& a_period,
                 const Vec3T<ShapePrecision<C>>& a_repeatLo,
                 const Vec3T<ShapePrecision<C>>& a_repeatHi);
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<BlurOp<ShapePrecision<C>>, C>
Blur(const C& a_operand, const ShapePrecision<C> a_blur, const ShapePrecision<C> a_alpha = ShapePrecision<C>(0.5));
 
template <class C>
[[nodiscard]] EBGEOMETRY_HOST
inline Node<MollifyOp<ShapePrecision<C>>, C>
Mollify(const C& a_operand, const ShapePrecision<C> a_dist, const size_t a_samples);
 
} // namespace Fused
 
} // namespace EBGeometry
 
#include "EBGeometry_FusedImplem.hpp"
 
// The ray forms of the fused operations and the BVH unions (and, through it, of the offsets): firstHit() is a
// template, and a file that included this header without them would march a composition that another file
// traces exactly, two definitions of one function. Included last, since that header includes this one.
#include "EBGeometry_CSGRayForms.hpp"
 
#endif