EBGeometry_SmoothRayForms.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_SMOOTHRAYFORMS_HPP
#define EBGEOMETRY_SMOOTHRAYFORMS_HPP
 
// Std includes
#include <cmath>
#include <cstdint>
#include <type_traits>
#include <utility>
 
// Our includes
#include "EBGeometry_CSGRayForms.hpp"
#include "EBGeometry_Fused.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_OffsetRayForms.hpp"
#include "EBGeometry_Ray.hpp"
#include "EBGeometry_RaySpans.hpp"
 
namespace EBGeometry {
 
namespace RayDetail {
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
blendDepth(const T a_smooth, const bool a_exponential) noexcept
{
  return a_exponential ? a_smooth * T(0.6931471805599453) : T(0.25) * a_smooth;
}
 
template <class S>
struct MarchedView
{
  using T = ShapePrecision<S>;
 
  static constexpr DistanceQuality distanceQuality = distanceQualityOf<S>;
 
  S m_shape;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    return m_shape.signedDistance(a_point);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return m_shape.computeBoundingVolume();
  }
};
 
template <class S>
struct IsComposition<MarchedView<S>> : std::bool_constant<CalledOutOfLine<S>>
{
};
 
template <class T, class... Cs>
struct BlendZoneView
{
  using Node = Fused::Node<Fused::SmoothCombineOp<T, true>, Cs...>;
 
  static constexpr DistanceQuality distanceQuality = distanceQualityOf<Node>;
 
  Node m_node;
 
  BoundingVolumes::AABBT<T> m_box;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    return this->value(a_point, std::index_sequence_for<Cs...>{});
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return m_box;
  }
 
private:
  template <size_t... Is>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  value(const Vec3T<T>& a_point, std::index_sequence<Is...>) const noexcept
  {
    const T smooth = m_node.layout().m_op.m_smooth;
 
    CSGOps::Pair<T> pair = CSGOps::pairInit<T>(false);
 
    (CSGOps::pairInsert(pair, RayDetail::fusedOperand<Is>(m_node).signedDistance(a_point), false), ...);
 
    const T blended = CSGOps::pairBlend(pair, smooth, false, false);
 
    return Math::max(blended, pair.m_second - (smooth + RayDetail::blendDepth(smooth, false)));
  }
};
 
template <class T, class... Cs>
struct IsComposition<BlendZoneView<T, Cs...>> : std::true_type
{
};
 
template <class T, class Op, class... Cs, size_t... Is>
EBGEOMETRY_HOST_DEVICE
inline void
pairOverlaps(const Fused::Node<Op, Cs...>& a_node,
             const T                       a_grow,
             BoundingVolumes::AABBT<T>&    a_box,
             bool&                         a_any,
             std::index_sequence<Is...>) noexcept
{
  const BoundingVolumes::AABBT<T> boxes[] = {RayDetail::fusedOperand<Is>(a_node).computeBoundingVolume()...};
  const Vec3T<T>                  grow    = a_grow * Vec3T<T>::ones();
 
  for (size_t i = 0; i < sizeof...(Is); i++) {
    for (size_t j = i + 1; j < sizeof...(Is); j++) {
      const Vec3T<T> lo = max(boxes[i].getLowCorner(), boxes[j].getLowCorner()) - grow;
      const Vec3T<T> hi = min(boxes[i].getHighCorner(), boxes[j].getHighCorner()) + grow;
 
      if (lo[0] <= hi[0] && lo[1] <= hi[1] && lo[2] <= hi[2]) {
        a_box = a_any ? BoundingVolumes::AABBT<T>(min(a_box.getLowCorner(), lo), max(a_box.getHighCorner(), hi))
                      : BoundingVolumes::AABBT<T>(lo, hi);
        a_any = true;
      }
    }
  }
}
 
} // namespace RayDetail
 
template <class T, bool IsUnion, class... Cs>
struct RayForm<Fused::Node<Fused::SmoothCombineOp<T, IsUnion>, Cs...>>
{
  using Node = Fused::Node<Fused::SmoothCombineOp<T, IsUnion>, Cs...>;
 
  [[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);
      };
 
    // A marched crossing's normal is the gradient there; marching again around it would find no better one.
    const auto resolve = [&a_shape, &a_ray, &a_options](const T a_t, const uint32_t) noexcept {
      const auto box = a_shape.computeBoundingVolume();
      const T    h   = a_options.m_tolerance * RayDetail::scaleOf(box.getLowCorner(), box.getHighCorner(), a_options);
 
      HitT<T> hit;
 
      hit.m_normal = RayDetail::gradientNormal(RayDetail::MarchedView<Node>{a_shape}, a_ray.at(a_t), h);
 
      return hit;
    };
 
    HitT<T> hit = RayDetail::traceSpans(trace, resolve, a_ray, a_tMin, a_tMax);
 
    hit.m_primitive = 0U;
 
    return hit;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline CrossingT<T>
  firstCrossing(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);
      };
 
    return RayDetail::traceCrossing(trace, a_ray, a_tMin, a_tMax);
  }
 
  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
  {
    RaySpansT<T> scratch;
    RaySpansT<T> operand;
 
    if constexpr (IsUnion) {
      const Fused::SmoothCombineOp<T, IsUnion>& op = a_shape.layout().m_op;
 
      if (op.m_exponential) {
        // An exponential blend differs from the plain union everywhere, if only slightly: marched whole.
        RayDetail::shapeSpans(
          RayDetail::MarchedView<Node>{a_shape}, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
      }
      else {
        // The plain union, exact, united with where the smooth union is inside and may differ from it, marched
        // up to what the plain union is known up to (see RayDetail::BlendZoneView).
        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);
 
        // The blend zone view, marched where two operands' boxes grown by w overlap, up to what the plain union is
        // known up to (see RayDetail::BlendZoneView).
        const T level = op.m_smooth + RayDetail::blendDepth(op.m_smooth, false);
 
        BoundingVolumes::AABBT<T> zone(Vec3T<T>::zeros(), Vec3T<T>::zeros());
 
        bool any = false;
 
        RayDetail::pairOverlaps(a_shape, level, zone, any, std::index_sequence_for<Cs...>{});
 
        // With no overlap, the smooth union is the plain union.
        if (any) {
          RayDetail::shapeSpans(RayDetail::BlendZoneView<T, Cs...>{a_shape, zone},
                                a_ray,
                                a_tMin,
                                a_spans.m_horizon,
                                a_options,
                                a_firstOnly,
                                operand);
          RayDetail::uniteInto(a_spans, operand);
 
          if (a_firstOnly) {
            RayDetail::keepFirst(a_spans);
          }
        }
      }
    }
    else {
      // Marched within the plain intersection, as a fused intersection traces it: the traced operands first, the
      // marched ones only within them. All its spans: where the smooth intersection is outside in one, it may be
      // inside in the next.
      constexpr size_t first = RayDetail::FirstTraced<T, Cs...>::value;
 
      RaySpansT<T> bracket;
 
      RayDetail::fusedOperandSpans<first>(a_shape, a_ray, a_tMin, a_tMax, a_options, false, bracket);
      RayDetail::fusedIntersect<0, false, first>(a_shape, a_ray, a_tMin, a_options, false, bracket, operand, scratch);
      RayDetail::fusedIntersect<0, true, first>(a_shape, a_ray, a_tMin, a_options, false, bracket, operand, scratch);
      RayDetail::spansWithin(
        RayDetail::MarchedView<Node>{a_shape}, a_ray, bracket, a_options, a_firstOnly, scratch, a_spans);
    }
 
    RayDetail::nameLeaf(0U, a_spans);
  }
};
 
} // namespace EBGeometry
 
#endif