EBGeometry_CSGRayForms.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_CSGRAYFORMS_HPP
#define EBGEOMETRY_CSGRAYFORMS_HPP
 
// Std includes
#include <cstddef>
#include <cstdint>
#include <type_traits>
#include <utility>
 
// Our includes
#include "EBGeometry_AnalyticRayForms.hpp"
#include "EBGeometry_BVHUnion.hpp"
#include "EBGeometry_CSGOps.hpp"
#include "EBGeometry_Fused.hpp"
#include "EBGeometry_MeshRayForms.hpp"
#include "EBGeometry_Ray.hpp"
 
namespace EBGeometry {
 
namespace Fused {
 
namespace Detail {
 
template <class Op>
struct TracedExactly : std::false_type
{
};
 
template <bool IsUnion>
struct TracedExactly<CombineOp<IsUnion>> : std::true_type
{
};
 
template <>
struct TracedExactly<ComplementOp> : std::true_type
{
};
 
template <class T>
struct TracedExactly<AffineOp<T>> : std::true_type
{
};
 
} // namespace Detail
 
} // namespace Fused
 
namespace RayDetail {
 
template <class S>
struct IsTracedFused : std::false_type
{
};
 
template <class Op, class... Cs>
struct IsTracedFused<Fused::Node<Op, Cs...>> : Fused::Detail::TracedExactly<Op>
{
};
 
template <class S>
struct LeafCount
{
  static constexpr uint32_t value = 1;
};
 
template <class Op, class... Cs>
struct LeafCount<Fused::Node<Op, Cs...>>
{
  static constexpr uint32_t value = Fused::Detail::TracedExactly<Op>::value ? (LeafCount<Cs>::value + ...) : 1U;
};
 
template <size_t I, class... Cs>
struct FirstLeaf;
 
template <class C, class... Cs>
struct FirstLeaf<0, C, Cs...>
{
  static constexpr uint32_t value = 0;
};
 
template <size_t I, class C, class... Cs>
struct FirstLeaf<I, C, Cs...>
{
  static constexpr uint32_t value = LeafCount<C>::value + FirstLeaf<I - 1, Cs...>::value;
};
 
template <class C, class = void>
struct DeferredOperand : std::false_type
{
};
 
template <class C, class T>
struct MarchedOperand
{
  static constexpr int value = (HasRayFirstHit<C, T>::value && !DeferredOperand<C>::value) ? 0 : 1;
};
 
template <class D, class T>
struct MarchedOperand<Fused::Node<Fused::ComplementOp, D>, T>
{
  static constexpr int value = (HasRayFirstHit<D, T>::value && !DeferredOperand<D>::value) ? 0 : 2;
};
 
template <class T, class... Cs>
struct FirstTraced
{
  static constexpr bool found = false;
 
  static constexpr size_t value = 0;
};
 
template <class T, class C, class... Cs>
struct FirstTraced<T, C, Cs...>
{
  static constexpr bool found = MarchedOperand<C, T>::value == 0 || FirstTraced<T, Cs...>::found;
 
  static constexpr size_t value =
    MarchedOperand<C, T>::value == 0 ? 0 : (FirstTraced<T, Cs...>::found ? 1 + FirstTraced<T, Cs...>::value : 0);
};
 
template <size_t I, class Op, class... Cs>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline auto
fusedOperand(const Fused::Node<Op, Cs...>& a_node) noexcept
{
  using C = typename Fused::Detail::TypeAt<I, Cs...>::type;
 
  return CSGDetail::fromStored<C>(Fused::Detail::get<I>(a_node.layout().m_operands), a_node.location());
}
 
template <size_t I, class T, class Op, class... Cs>
EBGEOMETRY_HOST_DEVICE
inline void
fusedOperandSpans(const Fused::Node<Op, Cs...>& a_node,
                  const RayT<T>&                a_ray,
                  const T                       a_t0,
                  const T                       a_t1,
                  const RayOptions<T>&          a_options,
                  const bool                    a_firstOnly,
                  RaySpansT<T>&                 a_spans) noexcept
{
  using C = typename Fused::Detail::TypeAt<I, Cs...>::type;
 
  RayDetail::shapeSpans(RayDetail::fusedOperand<I>(a_node), a_ray, a_t0, a_t1, a_options, a_firstOnly, a_spans);
 
  if constexpr (IsTracedFused<C>::value) {
    RayDetail::offsetLeaves(FirstLeaf<I, Cs...>::value, a_spans);
  }
  else {
    RayDetail::nameLeaf(FirstLeaf<I, Cs...>::value, a_spans);
  }
}
 
template <size_t I, class T, class Op, class... Cs>
EBGEOMETRY_HOST_DEVICE
inline void
fusedUnite(const Fused::Node<Op, Cs...>& a_node,
           const RayT<T>&                a_ray,
           const T                       a_t0,
           const RayOptions<T>&          a_options,
           const bool                    a_firstOnly,
           RaySpansT<T>&                 a_spans,
           RaySpansT<T>&                 a_operand) noexcept
{
  if constexpr (I < sizeof...(Cs)) {
    RayDetail::fusedOperandSpans<I>(a_node, a_ray, a_t0, a_spans.m_horizon, a_options, a_firstOnly, a_operand);
    RayDetail::uniteInto(a_spans, a_operand);
 
    if (a_firstOnly) {
      RayDetail::keepFirst(a_spans);
    }
 
    RayDetail::fusedUnite<I + 1>(a_node, a_ray, a_t0, a_options, a_firstOnly, a_spans, a_operand);
  }
}
 
template <size_t I, bool Marched, size_t Skip, class T, class Op, class... Cs>
EBGEOMETRY_HOST_DEVICE
inline void
fusedIntersect(const Fused::Node<Op, Cs...>& a_node,
               const RayT<T>&                a_ray,
               const T                       a_t0,
               const RayOptions<T>&          a_options,
               const bool                    a_firstOnly,
               RaySpansT<T>&                 a_spans,
               RaySpansT<T>&                 a_operand,
               RaySpansT<T>&                 a_scratch) noexcept
{
  if constexpr (I < sizeof...(Cs)) {
    using C = typename Fused::Detail::TypeAt<I, Cs...>::type;
 
    constexpr int kind = MarchedOperand<C, T>::value;
 
    if constexpr (I != Skip && (kind != 0) == Marched) {
      if constexpr (kind == 1) {
        RayDetail::spansWithin(
          RayDetail::fusedOperand<I>(a_node), a_ray, a_spans, a_options, a_firstOnly, a_scratch, a_operand);
        RayDetail::nameLeaf(FirstLeaf<I, Cs...>::value, a_operand);
        RayDetail::intersectInto(a_spans, a_operand, a_scratch);
      }
      else if constexpr (kind == 2) {
        // All the complemented operand's spans within the list: its complement there is the same. The trace
        // uses a_scratch as its scratch, which the complement then overwrites.
        RayDetail::spansWithin(RayDetail::fusedOperand<0>(RayDetail::fusedOperand<I>(a_node)),
                               a_ray,
                               a_spans,
                               a_options,
                               false,
                               a_scratch,
                               a_operand);
        RayDetail::nameLeaf(FirstLeaf<I, Cs...>::value, a_operand);
        RayDetail::complement(a_operand, a_t0, a_scratch);
        RayDetail::intersectInto(a_spans, a_scratch, a_operand);
      }
      else {
        RayDetail::fusedOperandSpans<I>(a_node, a_ray, a_t0, a_spans.m_horizon, a_options, false, a_operand);
        RayDetail::intersectInto(a_spans, a_operand, a_scratch);
      }
    }
 
    RayDetail::fusedIntersect<I + 1, Marched, Skip>(
      a_node, a_ray, a_t0, a_options, a_firstOnly, a_spans, a_operand, a_scratch);
  }
}
 
template <size_t I, class T, class Op, class... Cs>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline HitT<T>
fusedOperandLeafHit(const Fused::Node<Op, Cs...>& a_node,
                    const RayT<T>&                a_ray,
                    const T                       a_t,
                    const uint32_t                a_leaf,
                    const RayOptions<T>&          a_options) noexcept
{
  using C = typename Fused::Detail::TypeAt<I, Cs...>::type;
 
  if constexpr (IsTracedFused<C>::value) {
    return RayForm<C>::leafHit(
      RayDetail::fusedOperand<I>(a_node), a_ray, a_t, a_leaf - FirstLeaf<I, Cs...>::value, a_options);
  }
  else {
    return RayDetail::leafHit(RayDetail::fusedOperand<I>(a_node), a_ray, a_t, a_options);
  }
}
 
template <class T, class Op, class... Cs, size_t... Is>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline HitT<T>
fusedLeafHit(const Fused::Node<Op, Cs...>& a_node,
             const RayT<T>&                a_ray,
             const T                       a_t,
             const uint32_t                a_leaf,
             const RayOptions<T>&          a_options,
             std::index_sequence<Is...>) noexcept
{
  HitT<T> hit;
 
  const auto holds = [a_leaf](const uint32_t a_first, const uint32_t a_count) noexcept {
    return a_leaf >= a_first && a_leaf < a_first + a_count;
  };
 
  (void)((holds(FirstLeaf<Is, Cs...>::value, LeafCount<typename Fused::Detail::TypeAt<Is, Cs...>::type>::value)
            ? (hit = RayDetail::fusedOperandLeafHit<Is>(a_node, a_ray, a_t, a_leaf, a_options), true)
            : false) ||
         ...);
 
  return hit;
}
 
} // namespace RayDetail
 
template <class Op, class... Cs>
struct RayForm<Fused::Node<Op, Cs...>, std::enable_if_t<Fused::Detail::TracedExactly<Op>::value>>
{
  using Node = Fused::Node<Op, Cs...>;
 
  using T = typename Node::T;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(const Node& a_shape, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept
  {
    const auto trace =
      [&a_shape, &a_ray, &a_options](const T a_t0, const T a_t1, const bool a_all, RaySpansT<T>& a_spans) noexcept {
        spans(a_shape, a_ray, a_t0, a_t1, a_options, !a_all, a_spans);
      };
 
    const auto resolve = [&a_shape, &a_ray, &a_options](const T a_t, const uint32_t a_leaf) noexcept {
      return leafHit(a_shape, a_ray, a_t, a_leaf, a_options);
    };
 
    HitT<T> hit = RayDetail::traceSpans(trace, resolve, a_ray, a_tMin, a_tMax);
 
    hit.m_primitive = 0U;
 
    return hit;
  }
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const Node&          a_shape,
        const RayT<T>&       a_ray,
        T                    a_tMin,
        T                    a_tMax,
        const RayOptions<T>& a_options,
        bool                 a_firstOnly,
        RaySpansT<T>&        a_spans) noexcept
  {
    if constexpr (std::is_same_v<Op, Fused::UnionOp>) {
      RaySpansT<T> operand;
 
      RayDetail::fusedOperandSpans<0>(a_shape, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
 
      if (a_firstOnly) {
        RayDetail::keepFirst(a_spans);
      }
 
      RayDetail::fusedUnite<1>(a_shape, a_ray, a_tMin, a_options, a_firstOnly, a_spans, operand);
    }
    else if constexpr (std::is_same_v<Op, Fused::IntersectionOp>) {
      // An intersection needs all of each traced operand's spans: where the first one ends, the ray may still
      // be inside the others. The marched ones go last, marched only where the others leave the intersection
      // open.
      constexpr size_t first = RayDetail::FirstTraced<T, Cs...>::value;
 
      RaySpansT<T> operand;
      RaySpansT<T> scratch;
 
      RayDetail::fusedOperandSpans<first>(a_shape, a_ray, a_tMin, a_tMax, a_options, false, a_spans);
      RayDetail::fusedIntersect<0, false, first>(
        a_shape, a_ray, a_tMin, a_options, a_firstOnly, a_spans, operand, scratch);
      RayDetail::fusedIntersect<0, true, first>(
        a_shape, a_ray, a_tMin, a_options, a_firstOnly, a_spans, operand, scratch);
    }
    else if constexpr (std::is_same_v<Op, Fused::ComplementOp>) {
      RaySpansT<T> operand;
 
      RayDetail::fusedOperandSpans<0>(a_shape, a_ray, a_tMin, a_tMax, a_options, false, operand);
      RayDetail::complement(operand, a_tMin, a_spans);
    }
    else {
      const CSGOps::AffineCoefficients<T> map    = CSGOps::affineCoefficients(a_shape.layout().m_op.m_map);
      const RayDetail::MappedRay<T>       mapped = RayDetail::mapRay(map.m_values, a_ray, a_tMin, a_tMax);
 
      RayDetail::fusedOperandSpans<0>(a_shape, mapped.m_ray, mapped.m_t0, mapped.m_t1, a_options, a_firstOnly, a_spans);
      RayDetail::mapBack(mapped, a_tMin, a_tMax, a_spans);
    }
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  leafHit(const Node& a_shape, const RayT<T>& a_ray, T a_t, uint32_t a_leaf, const RayOptions<T>& a_options) noexcept
  {
    if constexpr (std::is_same_v<Op, Fused::UnionOp> || std::is_same_v<Op, Fused::IntersectionOp> ||
                  std::is_same_v<Op, Fused::ComplementOp>) {
      return RayDetail::fusedLeafHit(a_shape, a_ray, a_t, a_leaf, a_options, std::index_sequence_for<Cs...>{});
    }
    else {
      const CSGOps::AffineCoefficients<T> map    = CSGOps::affineCoefficients(a_shape.layout().m_op.m_map);
      const RayDetail::MappedRay<T>       mapped = RayDetail::mapRay(map.m_values, a_ray, a_t, a_t);
 
      HitT<T> hit = RayDetail::fusedLeafHit(
        a_shape, mapped.m_ray, mapped.m_t0, a_leaf, a_options, std::index_sequence_for<Cs...>{});
 
      hit.m_normal = RayDetail::mapNormalBack(map.m_values, hit.m_normal);
 
      return hit;
    }
  }
};
 
template <class T, class P, size_t K>
struct RayForm<BVHUnion<T, P, K>>
{
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(
    const BVHUnion<T, P, K>& a_shape, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept
  {
    const auto trace =
      [&a_shape, &a_ray, &a_options](const T a_t0, const T a_t1, const bool a_all, RaySpansT<T>& a_spans) noexcept {
        spans(a_shape, a_ray, a_t0, a_t1, a_options, !a_all, a_spans);
      };
 
    const auto resolve = [&a_shape, &a_ray, &a_options](const T a_t, const uint32_t a_leaf) noexcept {
      return leafHit(a_shape, a_ray, a_t, a_leaf, a_options);
    };
 
    return RayDetail::traceSpans(trace, resolve, a_ray, a_tMin, a_tMax);
  }
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const BVHUnion<T, P, K>& a_shape,
        const RayT<T>&           a_ray,
        T                        a_tMin,
        T                        a_tMax,
        const RayOptions<T>&     a_options,
        bool                     a_firstOnly,
        RaySpansT<T>&            a_spans) noexcept
  {
    using Root = typename BVHUnion<T, P, K>::Root;
 
    a_spans.reset(a_tMax);
 
    RaySpansT<T> operand;
 
    const PoolLocation location   = a_shape.location();
    const Root         bvh        = a_shape.getBVH();
    const auto         unbounded  = a_shape.getUnboundedPrimitives();
    const auto         primitives = bvh.getPrimitives();
 
    // Each primitive is traced only up to what the union so far is known up to, beyond which it can change
    // nothing known.
    const auto add = [&a_ray, a_tMin, &a_options, a_firstOnly, &location, &operand](
                       RaySpansT<T>& a_union, const CSGDetail::Stored<P>& a_primitive, const uint32_t a_leaf) noexcept {
      RayDetail::shapeSpans(CSGDetail::fromStored<P>(a_primitive, location),
                            a_ray,
                            a_tMin,
                            a_union.m_horizon,
                            a_options,
                            a_firstOnly,
                            operand);
      RayDetail::nameLeaf(a_leaf, operand);
      RayDetail::uniteInto(a_union, operand);
 
      if (a_firstOnly) {
        RayDetail::keepFirst(a_union);
      }
    };
 
    for (uint32_t i = 0; i < unbounded.size(); i++) {
      add(a_spans, unbounded[i], i);
    }
 
    const uint32_t numUnbounded = unbounded.size();
 
    const auto leaf =
      [&primitives, &add, numUnbounded](RaySpansT<T>& a_union, const size_t a_offset, const size_t a_count) noexcept {
        for (size_t i = a_offset; i < a_offset + a_count; i++) {
          add(a_union, primitives[static_cast<uint32_t>(i)], numUnbounded + static_cast<uint32_t>(i));
        }
      };
 
    // A box the ray enters beyond what the union is known up to can change nothing known.
    const auto bound = [](const RaySpansT<T>& a_union) noexcept -> T { return a_union.m_horizon; };
 
    bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), a_spans, leaf, bound);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  leafHit(const BVHUnion<T, P, K>& a_shape,
          const RayT<T>&           a_ray,
          T                        a_t,
          uint32_t                 a_leaf,
          const RayOptions<T>&     a_options) noexcept
  {
    const PoolLocation location  = a_shape.location();
    const auto         unbounded = a_shape.getUnboundedPrimitives();
 
    const auto primitives = a_shape.getBVH().getPrimitives();
 
    const CSGDetail::Stored<P> stored =
      (a_leaf < unbounded.size()) ? unbounded[a_leaf] : primitives[a_leaf - unbounded.size()];
 
    // The hit names the member crossed, not what the member would name (a mesh's face): a union's parts are
    // its members, as a tape's are its leaves.
    HitT<T> hit = RayDetail::leafHit(CSGDetail::fromStored<P>(stored, location), a_ray, a_t, a_options);
 
    hit.m_primitive = a_leaf;
 
    return hit;
  }
};
 
} // namespace EBGeometry
 
// The ray forms of the offsets and annuli, which build on these: whoever has the fused operations' ray forms
// has theirs too. Included last, since that header includes this one.
#include "EBGeometry_OffsetRayForms.hpp"
 
#endif