EBGeometry_OffsetRayForms.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_OFFSETRAYFORMS_HPP
#define EBGEOMETRY_OFFSETRAYFORMS_HPP
 
// Std includes
#include <cmath>
#include <cstdint>
#include <type_traits>
 
// Our includes
#include "EBGeometry_AnalyticDistanceFunctions.hpp"
#include "EBGeometry_AnalyticRayForms.hpp"
#include "EBGeometry_BVHQuery.hpp"
#include "EBGeometry_CSGRayForms.hpp"
#include "EBGeometry_DCEL_Iterator.hpp"
#include "EBGeometry_Fused.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_MeshDistanceFunctions.hpp"
#include "EBGeometry_MeshRayForms.hpp"
#include "EBGeometry_Ray.hpp"
#include "EBGeometry_RaySpans.hpp"
 
namespace EBGeometry {
 
namespace RayDetail {
 
template <class T>
struct Piece
{
  T m_in = -Math::Limits<T>::infinity();
 
  T m_out = Math::Limits<T>::infinity();
 
  Vec3T<T> m_nIn = Vec3T<T>::zeros();
 
  Vec3T<T> m_nOut = Vec3T<T>::zeros();
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline bool
  hit() const noexcept
  {
    return m_in <= m_out;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece
  empty() noexcept
  {
    Piece p;
 
    p.m_in  = Math::Limits<T>::infinity();
    p.m_out = -Math::Limits<T>::infinity();
 
    return p;
  }
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  clip(const RayT<T>& a_ray, const Vec3T<T>& a_n, const T a_d) noexcept
  {
    const T denom = dot(a_n, a_ray.m_direction);
    const T num   = a_d - dot(a_n, a_ray.m_origin);
 
    if (denom == T(0)) {
      if (num < T(0)) {
        *this = Piece::empty();
      }
 
      return;
    }
 
    const T t = num / denom;
 
    if (denom < T(0)) {
      if (t > m_in) {
        m_in  = t;
        m_nIn = a_n;
      }
    }
    else if (t < m_out) {
      m_out  = t;
      m_nOut = a_n;
    }
  }
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  unite(const Piece& a_other) noexcept
  {
    if (!a_other.hit()) {
      return;
    }
 
    if (!this->hit()) {
      *this = a_other;
 
      return;
    }
 
    if (a_other.m_in < m_in) {
      m_in  = a_other.m_in;
      m_nIn = a_other.m_nIn;
    }
 
    if (a_other.m_out > m_out) {
      m_out  = a_other.m_out;
      m_nOut = a_other.m_nOut;
    }
  }
};
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Piece<T>
ballPiece(const RayT<T>& a_ray, const Vec3T<T>& a_center, const T a_radius) noexcept
{
  const Vec3T<T> oc   = a_ray.m_origin - a_center;
  const T        b    = dot(oc, a_ray.m_direction);
  const Vec3T<T> perp = oc - b * a_ray.m_direction;
  const T        h    = a_radius * a_radius - dot(perp, perp);
 
  if (!(a_radius > T(0)) || h < T(0)) {
    return Piece<T>::empty();
  }
 
  const T s = std::sqrt(h);
 
  Piece<T> p;
 
  p.m_in   = -b - s;
  p.m_out  = -b + s;
  p.m_nIn  = (a_ray.at(p.m_in) - a_center) / a_radius;
  p.m_nOut = (a_ray.at(p.m_out) - a_center) / a_radius;
 
  return p;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Piece<T>
tubePiece(const RayT<T>& a_ray, const Vec3T<T>& a_point, const Vec3T<T>& a_axis, const T a_radius) noexcept
{
  if (!(a_radius > T(0))) {
    return Piece<T>::empty();
  }
 
  const Vec3T<T> oa = a_ray.m_origin - a_point;
  const Vec3T<T> dp = a_ray.m_direction - dot(a_ray.m_direction, a_axis) * a_axis;
  const Vec3T<T> op = oa - dot(oa, a_axis) * a_axis;
 
  const T qa = dot(dp, dp);
  const T qb = dot(dp, op);
 
  Piece<T> p;
 
  if (qa > T(0)) {
    const Vec3T<T> perp = op - (qb / qa) * dp;
    const T        h    = qa * (a_radius * a_radius - dot(perp, perp));
 
    if (h < T(0)) {
      return Piece<T>::empty();
    }
 
    const T s = std::sqrt(h);
 
    p.m_in   = (-qb - s) / qa;
    p.m_out  = (-qb + s) / qa;
    p.m_nIn  = (op + p.m_in * dp) / a_radius;
    p.m_nOut = (op + p.m_out * dp) / a_radius;
  }
  else if (dot(op, op) > a_radius * a_radius) {
    return Piece<T>::empty();
  }
 
  return p;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline Piece<T>
capsulePiece(const RayT<T>& a_ray, const Vec3T<T>& a_a, const Vec3T<T>& a_b, const T a_radius) noexcept
{
  Piece<T> p = RayDetail::ballPiece(a_ray, a_a, a_radius);
 
  p.unite(RayDetail::ballPiece(a_ray, a_b, a_radius));
 
  const Vec3T<T> axis   = a_b - a_a;
  const T        length = axis.length();
 
  if (length > T(0)) {
    const Vec3T<T> u = axis / length;
 
    Piece<T> tube = RayDetail::tubePiece(a_ray, a_a, u, a_radius);
 
    if (tube.hit()) {
      tube.clip(a_ray, -u, -dot(u, a_a));
      tube.clip(a_ray, u, dot(u, a_b));
 
      p.unite(tube);
    }
  }
 
  return p;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Piece<T>
boxPiece(const RayT<T>& a_ray, const Vec3T<T>& a_lo, const Vec3T<T>& a_hi) noexcept
{
  Piece<T> p;
 
  for (size_t dir = 0; dir < 3; dir++) {
    if (!(a_lo[dir] <= a_hi[dir])) {
      return Piece<T>::empty();
    }
 
    p.clip(a_ray, -Vec3T<T>::unit(dir), -a_lo[dir]);
    p.clip(a_ray, Vec3T<T>::unit(dir), a_hi[dir]);
  }
 
  return p;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline Piece<T>
roundedBoxPiece(const RayT<T>& a_ray, const Vec3T<T>& a_lo, const Vec3T<T>& a_hi, const T a_radius) noexcept
{
  // It lies in the box grown by the radius along every axis.
  const Vec3T<T> grow = a_radius * Vec3T<T>::ones();
 
  if (!RayDetail::boxPiece(a_ray, a_lo - grow, a_hi + grow).hit()) {
    return Piece<T>::empty();
  }
 
  // The three boxes grown along one axis each. Where the ray enters their union through a face across the axis
  // its box grows along, that face is the rounded box's, and the entry is the rounded box's (it is convex); so
  // for the exit. Near an edge it is a capsule's.
  Piece<T> p = Piece<T>::empty();
 
  size_t dirIn  = 0;
  size_t dirOut = 0;
 
  for (size_t dir = 0; dir < 3; dir++) {
    const Vec3T<T> along = a_radius * Vec3T<T>::unit(dir);
    const Piece<T> box   = RayDetail::boxPiece(a_ray, a_lo - along, a_hi + along);
 
    if (box.hit()) {
      if (!p.hit() || box.m_in < p.m_in) {
        dirIn = dir;
      }
 
      if (!p.hit() || box.m_out > p.m_out) {
        dirOut = dir;
      }
 
      p.unite(box);
    }
  }
 
  if (p.hit() && p.m_nIn[dirIn] != T(0) && p.m_nOut[dirOut] != T(0)) {
    return p;
  }
 
  for (size_t dir = 0; dir < 3; dir++) {
    const size_t d1 = (dir + 1) % 3;
    const size_t d2 = (dir + 2) % 3;
 
    for (int corner = 0; corner < 4; corner++) {
      Vec3T<T> a = a_lo;
 
      a[d1] = ((corner & 1) != 0) ? a_hi[d1] : a_lo[d1];
      a[d2] = ((corner & 2) != 0) ? a_hi[d2] : a_lo[d2];
 
      Vec3T<T> b = a;
 
      b[dir] = a_hi[dir];
 
      p.unite(RayDetail::capsulePiece(a_ray, a, b, a_radius));
    }
  }
 
  return p;
}
 
template <class T, class Edges>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline T
twiceArea(const Vec3T<T>& a_n, const Edges& a_edges) noexcept
{
  T        area  = T(0);
  bool     start = true;
  Vec3T<T> first = Vec3T<T>::zeros();
 
  a_edges([&](const Vec3T<T>& a_a, const Vec3T<T>& a_b) noexcept {
    if (start) {
      first = a_a;
      start = false;
    }
 
    area += dot(cross(a_a - first, a_b - a_a), a_n);
  });
 
  return area;
}
 
template <class T, class Edges>
EBGEOMETRY_HOST_DEVICE
inline bool
offsetFaceEntry(const RayT<T>&  a_ray,
                const Vec3T<T>& a_n,
                const Vec3T<T>& a_p0,
                const Edges&    a_edges,
                const T         a_radius,
                const T         a_tMin,
                const T         a_best,
                T&              a_t,
                Vec3T<T>&       a_normal,
                bool&           a_inSlab) noexcept
{
  const T dn = dot(a_n, a_ray.m_direction);
  const T h0 = dot(a_n, a_ray.m_origin - a_p0);
 
  T s0 = -Math::Limits<T>::infinity();
  T s1 = Math::Limits<T>::infinity();
 
  if (dn != T(0)) {
    const T a = (-a_radius - h0) / dn;
    const T b = (a_radius - h0) / dn;
 
    s0 = Math::min(a, b);
    s1 = Math::max(a, b);
  }
  else if (std::abs(h0) > a_radius) {
    a_inSlab = false;
 
    return false;
  }
 
  a_inSlab = s0 < a_tMin && s1 >= a_tMin;
 
  if (s1 < a_tMin || !(s0 < a_best)) {
    return false;
  }
 
  const T area = (dn != T(0)) ? RayDetail::twiceArea(a_n, a_edges) : T(0);
 
  if (area != T(0)) {
    // Edge values (edge length times the in-plane distance from the edge's line, positive inside) where the
    // ray's part in the slab begins and ends, projected onto the plane.
    const T        sign = (area > T(0)) ? T(1) : T(-1);
    const T        tIn  = Math::max(s0, a_tMin);
    const Vec3T<T> qIn  = a_ray.at(tIn) - dot(a_n, a_ray.at(tIn) - a_p0) * a_n;
    const Vec3T<T> qOut = a_ray.at(s1) - dot(a_n, a_ray.at(s1) - a_p0) * a_n;
 
    bool over      = true;
    bool separated = false;
 
    a_edges([&](const Vec3T<T>& a_a, const Vec3T<T>& a_b) noexcept {
      const Vec3T<T> e = a_b - a_a;
 
      const T in  = sign * dot(cross(e, qIn - a_a), a_n);
      const T out = sign * dot(cross(e, qOut - a_a), a_n);
      const T r   = a_radius * e.length();
 
      over      = over && in >= T(0);
      separated = separated || (in < -r && out < -r);
    });
 
    if (s0 >= a_tMin && over) {
      a_t      = s0;
      a_normal = ((dn < T(0)) ? T(1) : T(-1)) * a_n;
 
      return true;
    }
 
    if (separated) {
      return false;
    }
  }
 
  // The rim: each edge's capsule, if the ray passes within the radius of the edge's line.
  T    best  = a_best;
  bool found = false;
 
  a_edges([&](const Vec3T<T>& a_a, const Vec3T<T>& a_b) noexcept {
    const Vec3T<T> m     = cross(a_ray.m_direction, a_b - a_a);
    const T        m2    = dot(m, m);
    const T        along = dot(a_ray.m_origin - a_a, m);
 
    if (m2 > T(0) && along * along > a_radius * a_radius * m2) {
      return;
    }
 
    const Piece<T> capsule = RayDetail::capsulePiece(a_ray, a_a, a_b, a_radius);
 
    if (capsule.hit() && capsule.m_in >= a_tMin && capsule.m_in < best) {
      best     = capsule.m_in;
      a_t      = capsule.m_in;
      a_normal = capsule.m_nIn;
      found    = true;
    }
  });
 
  return found;
}
 
template <class T, class Edges>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline Piece<T>
offsetFacePiece(
  const RayT<T>& a_ray, const Vec3T<T>& a_n, const Vec3T<T>& a_p0, const Edges& a_edges, const T a_radius) noexcept
{
  Piece<T> p = Piece<T>::empty();
 
  const T area = RayDetail::twiceArea(a_n, a_edges);
 
  if (area != T(0)) {
    const T sign = (area > T(0)) ? T(1) : T(-1);
 
    Piece<T> prism;
 
    prism.clip(a_ray, a_n, dot(a_n, a_p0) + a_radius);
    prism.clip(a_ray, -a_n, -dot(a_n, a_p0) + a_radius);
 
    a_edges([&](const Vec3T<T>& a_a, const Vec3T<T>& a_b) noexcept {
      // Outward of the face, in its plane.
      const Vec3T<T> side = sign * cross(a_b - a_a, a_n);
      const T        len  = side.length();
 
      if (len > T(0)) {
        prism.clip(a_ray, side / len, dot(side / len, a_a));
      }
    });
 
    p = prism;
  }
 
  a_edges([&](const Vec3T<T>& a_a, const Vec3T<T>& a_b) noexcept {
    p.unite(RayDetail::capsulePiece(a_ray, a_a, a_b, a_radius));
  });
 
  return p;
}
 
template <class T>
struct GrownRayQuery
{
  using Scalar = T;
 
  using Key = T;
 
  Vec3T<T> m_origin;
 
  Vec3T<T> m_inverseDirection;
 
  T m_tMin;
 
  T m_grow;
 
  EBGEOMETRY_HOST_DEVICE
  inline GrownRayQuery(const RayT<T>& a_ray, const T a_tMin, const T a_grow) noexcept
    : m_origin(a_ray.m_origin),
      m_inverseDirection(T(1) / a_ray.m_direction[0], T(1) / a_ray.m_direction[1], T(1) / a_ray.m_direction[2]),
      m_tMin(a_tMin),
      m_grow(a_grow)
  {}
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  entry(const T (&a_lo)[3], const T (&a_hi)[3]) const noexcept
  {
    const T o[3]   = {m_origin[0], m_origin[1], m_origin[2]};
    const T inv[3] = {m_inverseDirection[0], m_inverseDirection[1], m_inverseDirection[2]};
    const T lo[3]  = {a_lo[0] - m_grow, a_lo[1] - m_grow, a_lo[2] - m_grow};
    const T hi[3]  = {a_hi[0] + m_grow, a_hi[1] + m_grow, a_hi[2] + m_grow};
 
    return BVH::QueryDetail::slabEntry<T>(lo, hi, o, inv, m_tMin);
  }
 
  template <size_t K>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline Array<T, K>
  keys(const BVH::ChildBoxes<T, K>& a_row) const noexcept
  {
    Array<T, K> keys;
 
    for (size_t k = 0; k < K; k++) {
      keys[k] = this->key(a_row, k);
    }
 
    return keys;
  }
 
  template <size_t K>
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  key(const BVH::ChildBoxes<T, K>& a_row, const size_t a_k) const noexcept
  {
    const T lo[3] = {a_row.m_lo[0][a_k], a_row.m_lo[1][a_k], a_row.m_lo[2][a_k]};
    const T hi[3] = {a_row.m_hi[0][a_k], a_row.m_hi[1][a_k], a_row.m_hi[2][a_k]};
 
    return this->entry(lo, hi);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  key(const BoundingVolumes::AABBT<T>& a_box) const noexcept
  {
    const Vec3T<T>& l = a_box.getLowCorner();
    const Vec3T<T>& h = a_box.getHighCorner();
 
    const T lo[3] = {l[0], l[1], l[2]};
    const T hi[3] = {h[0], h[1], h[2]};
 
    return this->entry(lo, hi);
  }
};
 
template <class T>
struct OffsetFaceUnion
{
  static constexpr uint32_t Capacity = 4;
 
  T m_in[Capacity];
 
  T m_out[Capacity];
 
  Vec3T<T> m_nOut[Capacity];
 
  uint32_t m_count = 0;
 
  T m_horizon = Math::Limits<T>::infinity();
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  insert(T a_in, T a_out, Vec3T<T> a_nOut) noexcept
  {
    if (a_in > m_horizon) {
      return;
    }
 
    // The spans it meets, [first, last), merged into one.
    uint32_t first = 0;
 
    while (first < m_count && m_out[first] < a_in) {
      first++;
    }
 
    uint32_t last = first;
 
    while (last < m_count && m_in[last] <= a_out) {
      a_in = Math::min(a_in, m_in[last]);
 
      if (m_out[last] > a_out) {
        a_out  = m_out[last];
        a_nOut = m_nOut[last];
      }
 
      last++;
    }
 
    if (last > first) {
      m_in[first]   = a_in;
      m_out[first]  = a_out;
      m_nOut[first] = a_nOut;
 
      const uint32_t removed = last - first - 1U;
 
      for (uint32_t k = first + 1U; k + removed < m_count; k++) {
        m_in[k]   = m_in[k + removed];
        m_out[k]  = m_out[k + removed];
        m_nOut[k] = m_nOut[k + removed];
      }
 
      m_count -= removed;
 
      return;
    }
 
    // A new span at position first.
    if (m_count == Capacity) {
      if (first == m_count) {
        m_horizon = a_in;
 
        return;
      }
 
      m_horizon = m_in[Capacity - 1U];
 
      m_count--;
    }
 
    for (uint32_t k = m_count; k > first; k--) {
      m_in[k]   = m_in[k - 1U];
      m_out[k]  = m_out[k - 1U];
      m_nOut[k] = m_nOut[k - 1U];
    }
 
    m_in[first]   = a_in;
    m_out[first]  = a_out;
    m_nOut[first] = a_nOut;
 
    m_count++;
  }
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  toSpans(const T a_t0, RaySpansT<T>& a_spans) const noexcept
  {
    a_spans.reset(m_horizon);
 
    for (uint32_t i = 0; i < m_count && m_in[i] <= m_horizon; i++) {
      const SpanEndT<T> a = (m_in[i] < a_t0) ? RayDetail::openAt(a_t0) : SpanEndT<T>{m_in[i], 0U, SpanEndT<T>::Exact};
      const SpanEndT<T> b =
        (m_out[i] >= m_horizon) ? RayDetail::openAt(m_horizon) : SpanEndT<T>{m_out[i], 0U, SpanEndT<T>::Exact};
 
      a_spans.add(a, b);
    }
  }
};
 
template <class S>
struct MeshFaces
{
};
 
template <class T, size_t K, size_t W>
struct MeshFaces<TriMeshSDF<T, K, W>>
{
  template <class State, class Face, class Bound>
  EBGEOMETRY_HOST_DEVICE
  static inline void
  walk(const TriMeshSDF<T, K, W>& a_mesh,
       const RayT<T>&             a_ray,
       const T                    a_tMin,
       const T                    a_grow,
       State&                     a_state,
       const Face&                a_face,
       const Bound&               a_bound) noexcept
  {
    const auto bvh    = a_mesh.getRoot();
    const auto groups = bvh.getPrimitives();
 
    const auto leaf = [&groups, &a_face](State& a_s, size_t a_offset, size_t a_count) noexcept {
      for (size_t i = a_offset; i < a_offset + a_count; i++) {
        const auto& group = groups[static_cast<uint32_t>(i)];
 
        for (size_t lane = 0; lane < W; lane++) {
          const Vec3T<T> v[3] = {group.getVertex(lane, 0), group.getVertex(lane, 1), group.getVertex(lane, 2)};
 
          const auto edges = [&v](const auto& a_f) noexcept {
            a_f(v[0], v[1]);
            a_f(v[1], v[2]);
            a_f(v[2], v[0]);
          };
 
          const auto distance2 = [&group, lane](const Vec3T<T>& a_point) noexcept {
            const T d = group.signedDistanceAt(static_cast<uint32_t>(lane), a_point);
 
            return d * d;
          };
 
          a_face(a_s, group.getNormal(lane), v[0], edges, distance2, group.getFaceId(lane));
        }
      }
    };
 
    bvh.traverse(GrownRayQuery<T>(a_ray, a_tMin, a_grow), a_state, leaf, a_bound);
  }
};
 
template <class T, size_t K>
struct MeshFaces<MeshSDF<T, K>>
{
  template <class State, class Face, class Bound>
  EBGEOMETRY_HOST_DEVICE
  static inline void
  walk(const MeshSDF<T, K>& a_mesh,
       const RayT<T>&       a_ray,
       const T              a_tMin,
       const T              a_grow,
       State&               a_state,
       const Face&          a_face,
       const Bound&         a_bound) noexcept
  {
    const auto bvh     = a_mesh.getRoot();
    const auto mesh    = a_mesh.getMesh();
    const auto faceIds = bvh.getPrimitives();
 
    const auto leaf = [&faceIds, &mesh, &a_face](State& a_s, size_t a_offset, size_t a_count) noexcept {
      for (size_t i = a_offset; i < a_offset + a_count; i++) {
        const uint32_t id = faceIds[static_cast<uint32_t>(i)];
 
        MeshFaces::visit(mesh, id, a_s, a_face);
      }
    };
 
    bvh.traverse(GrownRayQuery<T>(a_ray, a_tMin, a_grow), a_state, leaf, a_bound);
  }
 
  template <class Mesh, class State, class Face>
  EBGEOMETRY_HOST_DEVICE
  static inline void
  visit(const Mesh& a_mesh, const uint32_t a_id, State& a_state, const Face& a_face) noexcept
  {
    const auto& face = a_mesh.getFace(a_id);
 
    const auto edges = [&a_mesh, &face](const auto& a_f) noexcept {
      for (DCEL::EdgeIteratorT<T> iter(a_mesh, face); iter.ok(); ++iter) {
        const auto& edge = a_mesh.getEdge(iter());
 
        a_f(edge.getVertex(a_mesh).getPosition(), edge.getNextEdge(a_mesh).getVertex(a_mesh).getPosition());
      }
    };
 
    const auto distance2 = [&a_mesh, &face](const Vec3T<T>& a_point) noexcept {
      return face.unsignedDistance2(a_point, a_mesh);
    };
 
    a_face(a_state, face.getNormal(), face.getCentroid(), edges, distance2, a_id);
  }
};
 
template <class T>
struct MeshFaces<FlatMeshSDF<T>>
{
  template <class State, class Face, class Bound>
  EBGEOMETRY_HOST_DEVICE
  static inline void
  walk(const FlatMeshSDF<T>& a_mesh,
       const RayT<T>& /*a_ray*/,
       const T /*a_tMin*/,
       const T /*a_grow*/,
       State&      a_state,
       const Face& a_face,
       const Bound& /*a_bound*/) noexcept
  {
    const auto mesh = a_mesh.getMesh();
 
    for (uint32_t i = 0; i < mesh.numFaces(); i++) {
      MeshFaces<MeshSDF<T, 4>>::visit(mesh, i, a_state, a_face);
    }
  }
};
 
template <class T>
struct NearestOffsetFace
{
  T m_t;
 
  uint32_t m_face;
 
  Vec3T<T> m_normal;
 
  bool m_inside;
};
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline NearestOffsetFace<T>
nearestOffsetFace(const S& a_mesh, const RayT<T>& a_ray, const T a_tMin, const T a_tMax, const T a_radius) noexcept
{
  NearestOffsetFace<T> found{a_tMax, UINT32_MAX, Vec3T<T>::zeros(), false};
 
  const auto face = [&a_ray, a_tMin, a_radius](NearestOffsetFace<T>& a_state,
                                               const Vec3T<T>&       a_n,
                                               const Vec3T<T>&       a_p0,
                                               const auto&           a_edges,
                                               const auto&           a_distance2,
                                               const uint32_t        a_id) noexcept {
    if (a_state.m_inside) {
      return;
    }
 
    T        t      = T(0);
    Vec3T<T> normal = Vec3T<T>::zeros();
    bool     inSlab = false;
 
    const bool entered =
      RayDetail::offsetFaceEntry(a_ray, a_n, a_p0, a_edges, a_radius, a_tMin, a_state.m_t, t, normal, inSlab);
 
    // A start within the radius of the face is inside the grown surface.
    if (inSlab && a_distance2(a_ray.at(a_tMin)) <= a_radius * a_radius) {
      a_state.m_inside = true;
      a_state.m_t      = a_tMin;
 
      return;
    }
 
    if (entered) {
      a_state.m_t      = t;
      a_state.m_face   = a_id;
      a_state.m_normal = normal;
    }
  };
 
  const auto bound = [](const NearestOffsetFace<T>& a_state) noexcept -> T { return a_state.m_t; };
 
  MeshFaces<S>::walk(a_mesh, a_ray, a_tMin, a_radius, found, face, bound);
 
  return found;
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline void
offsetFaceUnion(const S&            a_mesh,
                const RayT<T>&      a_ray,
                const T             a_tMin,
                const T             a_tMax,
                const T             a_radius,
                OffsetFaceUnion<T>& a_union) noexcept
{
  a_union.m_count   = 0;
  a_union.m_horizon = a_tMax;
 
  const auto face = [&a_ray, a_tMin, a_radius](OffsetFaceUnion<T>& a_state,
                                               const Vec3T<T>&     a_n,
                                               const Vec3T<T>&     a_p0,
                                               const auto&         a_edges,
                                               const auto& /*a_distance2*/,
                                               const uint32_t /*a_id*/) noexcept {
    const Piece<T> p = RayDetail::offsetFacePiece(a_ray, a_n, a_p0, a_edges, a_radius);
 
    if (p.hit() && p.m_out >= a_tMin && p.m_in <= a_state.m_horizon) {
      a_state.insert(p.m_in, p.m_out, p.m_nOut);
    }
  };
 
  const auto bound = [](const OffsetFaceUnion<T>& a_state) noexcept -> T { return a_state.m_horizon; };
 
  MeshFaces<S>::walk(a_mesh, a_ray, a_tMin, a_radius, a_union, face, bound);
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline void
meshSpans(const S&             a_mesh,
          const RayT<T>&       a_ray,
          const T              a_tMin,
          const T              a_tMax,
          const RayOptions<T>& a_options,
          const bool           a_firstOnly,
          RaySpansT<T>&        a_spans) noexcept
{
  RayDetail::shapeSpans(a_mesh, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
}
 
template <class S>
struct IsMeshLevel : std::false_type
{
};
 
template <class T, size_t K, size_t W>
struct IsMeshLevel<TriMeshSDF<T, K, W>> : std::true_type
{
};
 
template <class T, size_t K>
struct IsMeshLevel<MeshSDF<T, K>> : std::true_type
{
};
 
template <class T>
struct IsMeshLevel<FlatMeshSDF<T>> : std::true_type
{
};
 
template <class T, class C>
struct IsMeshLevel<Fused::Node<Fused::AffineOp<T>, C>> : IsMeshLevel<C>
{
};
 
template <class T, bool IsAnnulus, class C>
struct DeferredOperand<Fused::Node<Fused::ShiftOp<T, IsAnnulus>, C>, std::enable_if_t<IsMeshLevel<C>::value>>
  : std::true_type
{
};
 
template <class S, class = void>
struct LevelForm
{
};
 
template <class T>
struct LevelForm<SphereSDF<T>>
{
  static constexpr bool Convex = true;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece<T>
  piece(const SphereSDF<T>& a_shape, const T a_level, const RayT<T>& a_ray) noexcept
  {
    return RayDetail::ballPiece(a_ray, a_shape.getCenter(), a_shape.getRadius() + a_level);
  }
};
 
template <class T>
struct LevelForm<CapsuleSDF<T>>
{
  static constexpr bool Convex = true;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece<T>
  piece(const CapsuleSDF<T>& a_shape, const T a_level, const RayT<T>& a_ray) noexcept
  {
    return RayDetail::capsulePiece(a_ray, a_shape.getCenter1(), a_shape.getCenter2(), a_shape.getRadius() + a_level);
  }
};
 
template <class T>
struct LevelForm<InfiniteCylinderSDF<T>>
{
  static constexpr bool Convex = true;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece<T>
  piece(const InfiniteCylinderSDF<T>& a_shape, const T a_level, const RayT<T>& a_ray) noexcept
  {
    return RayDetail::tubePiece(
      a_ray, a_shape.getCenter(), Vec3T<T>::unit(a_shape.getAxis()), a_shape.getRadius() + a_level);
  }
};
 
template <class T>
struct LevelForm<PlaneSDF<T>>
{
  static constexpr bool Convex = true;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece<T>
  piece(const PlaneSDF<T>& a_shape, const T a_level, const RayT<T>& a_ray) noexcept
  {
    Piece<T> p;
 
    p.clip(a_ray, a_shape.getNormal(), dot(a_shape.getNormal(), a_shape.getPoint()) + a_level);
 
    return p;
  }
};
 
template <class T>
struct LevelForm<BoxSDF<T>>
{
  static constexpr bool Convex = true;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline Piece<T>
  piece(const BoxSDF<T>& a_shape, const T a_level, const RayT<T>& a_ray) noexcept
  {
    const Vec3T<T>& lo = a_shape.getLowCorner();
    const Vec3T<T>& hi = a_shape.getHighCorner();
 
    if (a_level > T(0)) {
      return RayDetail::roundedBoxPiece(a_ray, lo, hi, a_level);
    }
 
    const Vec3T<T> shrink = -a_level * Vec3T<T>::ones();
 
    return RayDetail::boxPiece(a_ray, lo + shrink, hi - shrink);
  }
};
 
template <class S>
struct MeshLevelForm
{
  static constexpr bool Convex = false;
 
  static constexpr bool Mesh = true;
};
 
template <class T, size_t K, size_t W>
struct LevelForm<TriMeshSDF<T, K, W>> : MeshLevelForm<TriMeshSDF<T, K, W>>
{
};
 
template <class T, size_t K>
struct LevelForm<MeshSDF<T, K>> : MeshLevelForm<MeshSDF<T, K>>
{
};
 
template <class T>
struct LevelForm<FlatMeshSDF<T>> : MeshLevelForm<FlatMeshSDF<T>>
{
};
 
template <class S, class = void>
struct HasLevelForm : std::false_type
{
};
 
template <class S>
struct HasLevelForm<S, std::void_t<decltype(LevelForm<S>::Convex)>> : std::true_type
{
};
 
template <class T, class C>
struct LevelForm<Fused::Node<Fused::AffineOp<T>, C>, std::enable_if_t<HasLevelForm<C>::value>>
{
  static constexpr bool Convex = false;
 
  static constexpr bool Mesh = false;
};
 
template <class S>
struct LevelView
{
  using T = ShapePrecision<S>;
 
  static constexpr DistanceQuality distanceQuality = distanceQualityOf<S>;
 
  S m_shape;
 
  T m_level;
 
  bool m_annulus;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    const T f = m_shape.signedDistance(a_point);
 
    return (m_annulus ? std::abs(f) : f) - m_level;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    const BoundingVolumes::AABBT<T> box  = m_shape.computeBoundingVolume();
    const Vec3T<T>                  grow = Math::max(m_level, T(0)) * Vec3T<T>::ones();
 
    return BoundingVolumes::AABBT<T>(box.getLowCorner() - grow, box.getHighCorner() + grow);
  }
};
 
template <class S>
struct IsComposition<LevelView<S>> : std::bool_constant<CalledOutOfLine<S>>
{
};
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
inline int
convexLevelSpans(const LevelView<S>& a_view, const RayT<T>& a_ray, Piece<T> (&a_spans)[2]) noexcept
{
  if (a_view.m_annulus && !(a_view.m_level > T(0))) {
    return 0;
  }
 
  const Piece<T> outer = LevelForm<S>::piece(a_view.m_shape, a_view.m_level, a_ray);
 
  if (!outer.hit()) {
    return 0;
  }
 
  const Piece<T> inner =
    a_view.m_annulus ? LevelForm<S>::piece(a_view.m_shape, -a_view.m_level, a_ray) : Piece<T>::empty();
 
  if (!inner.hit()) {
    a_spans[0] = outer;
 
    return 1;
  }
 
  int n = 0;
 
  if (inner.m_in > outer.m_in) {
    a_spans[n++] = Piece<T>{outer.m_in, inner.m_in, outer.m_nIn, -inner.m_nIn};
  }
 
  if (outer.m_out > inner.m_out) {
    a_spans[n++] = Piece<T>{inner.m_out, outer.m_out, -inner.m_nOut, outer.m_nOut};
  }
 
  return n;
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline HitT<T>
levelHitAt(
  const RayT<T>& a_ray, const T a_t, const Vec3T<T>& a_normal, const bool a_entering, const uint32_t a_primitive)
{
  HitT<T> hit;
 
  hit.m_status    = HitStatus::Hit;
  hit.m_t         = a_t;
  hit.m_point     = a_ray.at(a_t);
  hit.m_normal    = a_normal;
  hit.m_entering  = a_entering;
  hit.m_exact     = true;
  hit.m_primitive = a_primitive;
 
  return hit;
}
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
levelNormal(const LevelView<S>& a_view, const Vec3T<T>& a_point, const T a_h) noexcept
{
  const Vec3T<T> n = RayDetail::gradientNormal(a_view.m_shape, a_point, a_h);
 
  return (a_view.m_annulus && RayDetail::valueAt(a_view.m_shape, a_point) < T(0)) ? -n : n;
}
 
// A level view's spans, out of line on a device: declared here for its ray form, defined after it.
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline void
levelSpans(const LevelView<S>&  a_view,
           const RayT<T>&       a_ray,
           T                    a_tMin,
           T                    a_tMax,
           const RayOptions<T>& a_options,
           bool                 a_firstOnly,
           RaySpansT<T>&        a_spans) noexcept;
 
} // namespace RayDetail
 
template <class S>
struct RayForm<RayDetail::LevelView<S>, std::enable_if_t<RayDetail::HasLevelForm<S>::value>>
{
  using View = RayDetail::LevelView<S>;
 
  using T = typename View::T;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(const View& a_view, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept
  {
    using Form = RayDetail::LevelForm<S>;
 
    if constexpr (Form::Convex) {
      RayDetail::Piece<T> spans[2];
 
      const int n = RayDetail::convexLevelSpans(a_view, a_ray, spans);
 
      for (int i = 0; i < n; i++) {
        const RayDetail::Piece<T>& s = spans[i];
 
        if (s.m_in >= a_tMin && s.m_in <= a_tMax) {
          return RayDetail::levelHitAt(a_ray, s.m_in, s.m_nIn, true, 0U);
        }
 
        if (s.m_in < a_tMin && s.m_out >= a_tMin && s.m_out <= a_tMax) {
          return RayDetail::levelHitAt(a_ray, s.m_out, s.m_nOut, false, 0U);
        }
      }
 
      return HitT<T>();
    }
    else if constexpr (Form::Mesh) {
      if (a_view.m_annulus || a_view.m_level == T(0)) {
        if (a_view.m_level == T(0) && !a_view.m_annulus) {
          return RayDetail::shapeHit(a_view.m_shape, a_ray, a_tMin, a_tMax, a_options);
        }
 
        return annulusHit(a_view, a_ray, a_tMin, a_tMax);
      }
 
      // An offset of a mesh: from its spans, the crossing's normal from the mesh (leaf 0) or its offset faces
      // (leaf 1).
      const auto trace =
        [&a_view, &a_ray, &a_options](const T a_t0, const T a_t1, const bool a_all, RaySpansT<T>& a_spans) noexcept {
          RayDetail::levelSpans(a_view, a_ray, a_t0, a_t1, a_options, !a_all, a_spans);
        };
 
      const auto resolve = [&a_view, &a_ray, &a_options](const T a_t, const uint32_t a_leaf) noexcept {
        if (a_leaf == 0U) {
          return RayDetail::leafHit(a_view.m_shape, a_ray, a_t, a_options);
        }
 
        return annulusLeafHit(View{a_view.m_shape, std::abs(a_view.m_level), true}, a_ray, a_t, a_options);
      };
 
      return RayDetail::traceSpans(trace, resolve, a_ray, a_tMin, a_tMax);
    }
    else {
      const auto& map = a_view.m_shape.layout().m_op;
 
      const CSGOps::AffineCoefficients<T> m      = CSGOps::affineCoefficients(map.m_map);
      const RayDetail::MappedRay<T>       mapped = RayDetail::mapRay(m.m_values, a_ray, a_tMin, a_tMax);
 
      const auto inner = mappedView(a_view);
 
      HitT<T> hit = EBGeometry::firstHit(inner, mapped.m_ray, mapped.m_t0, mapped.m_t1, a_options);
 
      if (hit.m_status != HitStatus::Miss) {
        hit.m_t = Math::clamp(hit.m_t / mapped.m_scale, a_tMin, a_tMax);
      }
 
      if (hit.hit()) {
        hit.m_point  = a_ray.at(hit.m_t);
        hit.m_normal = RayDetail::mapNormalBack(m.m_values, hit.m_normal);
      }
 
      return hit;
    }
  }
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const View&          a_view,
        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 Form = RayDetail::LevelForm<S>;
 
    if constexpr (Form::Convex) {
      RayDetail::Piece<T> spans[2];
 
      const int n = RayDetail::convexLevelSpans(a_view, a_ray, spans);
 
      a_spans.reset(a_tMax);
 
      for (int i = 0; i < n; i++) {
        a_spans.addClipped(spans[i].m_in, spans[i].m_out, a_tMin, a_tMax);
      }
    }
    else if constexpr (Form::Mesh) {
      const T level = a_view.m_level;
 
      if (!a_view.m_annulus && level == T(0)) {
        RayDetail::meshSpans(a_view.m_shape, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
 
        return;
      }
 
      if (a_view.m_annulus) {
        annulusSpans(a_view, a_ray, a_tMin, a_tMax, a_firstOnly, a_spans);
 
        return;
      }
 
      // An offset: the mesh (leaf 0) united with its faces offset by the level (leaf 1), or less those offset
      // by its opposite.
      const View   faces{a_view.m_shape, std::abs(level), true};
      RaySpansT<T> grown;
      RaySpansT<T> scratch;
 
      if (level > T(0)) {
        RayDetail::meshSpans(a_view.m_shape, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
        annulusSpans(faces, a_ray, a_tMin, a_spans.m_horizon, a_firstOnly, grown);
        RayDetail::nameLeaf(1U, grown);
        RayDetail::uniteInto(a_spans, grown);
 
        if (a_firstOnly) {
          RayDetail::keepFirst(a_spans);
        }
      }
      else {
        RayDetail::meshSpans(a_view.m_shape, a_ray, a_tMin, a_tMax, a_options, false, a_spans);
        annulusSpans(faces, a_ray, a_tMin, a_spans.m_horizon, false, grown);
        RayDetail::nameLeaf(1U, grown);
        RayDetail::complement(grown, a_tMin, scratch);
        RayDetail::intersectInto(a_spans, scratch, grown);
      }
    }
    else {
      const auto& map = a_view.m_shape.layout().m_op;
 
      const CSGOps::AffineCoefficients<T> m      = CSGOps::affineCoefficients(map.m_map);
      const RayDetail::MappedRay<T>       mapped = RayDetail::mapRay(m.m_values, a_ray, a_tMin, a_tMax);
 
      RayDetail::shapeSpans(
        mappedView(a_view), mapped.m_ray, mapped.m_t0, mapped.m_t1, a_options, a_firstOnly, a_spans);
      RayDetail::mapBack(mapped, a_tMin, a_tMax, a_spans);
    }
  }
 
private:
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline auto
  mappedView(const View& a_view) noexcept
  {
    const auto& op      = a_view.m_shape.layout().m_op;
    const auto  operand = RayDetail::fusedOperand<0>(a_view.m_shape);
 
    return RayDetail::LevelView<std::decay_t<decltype(operand)>>{
      operand, a_view.m_level / op.m_valueScale, a_view.m_annulus};
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  annulusHit(const View& a_view, const RayT<T>& a_ray, const T a_tMin, const T a_tMax) noexcept
  {
    if (!(a_view.m_level > T(0))) {
      return HitT<T>();
    }
 
    const RayDetail::NearestOffsetFace<T> found =
      RayDetail::nearestOffsetFace(a_view.m_shape, a_ray, a_tMin, a_tMax, a_view.m_level);
 
    if (!found.m_inside) {
      return (found.m_face == UINT32_MAX) ? HitT<T>()
                                          : RayDetail::levelHitAt(a_ray, found.m_t, found.m_normal, true, found.m_face);
    }
 
    RayDetail::OffsetFaceUnion<T> merged;
 
    RayDetail::offsetFaceUnion(a_view.m_shape, a_ray, a_tMin, a_tMax, a_view.m_level, merged);
 
    if (merged.m_count > 0U && merged.m_out[0] < merged.m_horizon) {
      return RayDetail::levelHitAt(a_ray, merged.m_out[0], merged.m_nOut[0], false, 0U);
    }
 
    HitT<T> hit;
 
    if (merged.m_horizon < a_tMax) {
      hit.m_status = HitStatus::Unresolved;
      hit.m_t      = merged.m_horizon;
    }
 
    return hit;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  annulusLeafHit(const View& a_view, const RayT<T>& a_ray, const T a_t, const RayOptions<T>& a_options) noexcept
  {
    const auto box   = a_view.computeBoundingVolume();
    const T    width = a_options.m_tolerance * RayDetail::scaleOf(box.getLowCorner(), box.getHighCorner(), a_options);
 
    HitT<T> hit = annulusHit(a_view, a_ray, a_t - width, a_t + width);
 
    if (!hit.hit()) {
      hit.m_normal = RayDetail::levelNormal(a_view, a_ray.at(a_t), width);
    }
 
    return hit;
  }
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  annulusSpans(const View&    a_view,
               const RayT<T>& a_ray,
               const T        a_tMin,
               const T        a_tMax,
               const bool     a_firstOnly,
               RaySpansT<T>&  a_spans) noexcept
  {
    a_spans.reset(a_tMax);
 
    if (!(a_view.m_level > T(0))) {
      return;
    }
 
    if (a_firstOnly) {
      const RayDetail::NearestOffsetFace<T> found =
        RayDetail::nearestOffsetFace(a_view.m_shape, a_ray, a_tMin, a_tMax, a_view.m_level);
 
      if (!found.m_inside) {
        if (found.m_face != UINT32_MAX) {
          a_spans.add(SpanEndT<T>{found.m_t, 0U, SpanEndT<T>::Exact}, RayDetail::openAt(found.m_t));
 
          a_spans.m_horizon = found.m_t;
        }
 
        return;
      }
    }
 
    RayDetail::OffsetFaceUnion<T> merged;
 
    RayDetail::offsetFaceUnion(a_view.m_shape, a_ray, a_tMin, a_tMax, a_view.m_level, merged);
 
    merged.toSpans(a_tMin, a_spans);
  }
};
 
namespace RayDetail {
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline void
levelSpans(const LevelView<S>&  a_view,
           const RayT<T>&       a_ray,
           const T              a_tMin,
           const T              a_tMax,
           const RayOptions<T>& a_options,
           const bool           a_firstOnly,
           RaySpansT<T>&        a_spans) noexcept
{
  RayForm<LevelView<S>>::spans(a_view, a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline void
levelSpansWithin(const LevelView<S>&  a_view,
                 const RayT<T>&       a_ray,
                 const RaySpansT<T>&  a_within,
                 const T              a_scale,
                 const RayOptions<T>& a_options,
                 const bool           a_firstOnly,
                 RaySpansT<T>&        a_spans) noexcept
{
  const T horizon = a_within.m_horizon * a_scale;
 
  a_spans.reset(horizon);
 
  RaySpansT<T> piece;
 
  for (uint32_t i = 0; i < a_within.m_count && a_within.m_spans[i].m_in.m_t * a_scale <= horizon; i++) {
    const T t0 = a_within.m_spans[i].m_in.m_t * a_scale;
    const T t1 = Math::min(a_within.m_spans[i].m_out.m_t * a_scale, horizon);
 
    RayDetail::levelSpans(a_view, a_ray, t0, t1, a_options, a_firstOnly, piece);
 
    for (uint32_t k = 0; k < piece.m_count && piece.m_spans[k].m_in.m_t <= piece.m_horizon; k++) {
      const typename RaySpansT<T>::Span span = RayDetail::cut(piece.m_spans[k], piece.m_horizon);
 
      a_spans.add(span.m_in, span.m_out);
    }
 
    if (piece.m_horizon < t1 || (a_firstOnly && piece.m_count > 0U)) {
      a_spans.m_horizon = Math::min(a_spans.m_horizon, piece.m_horizon);
    }
 
    // Cut back to this span or before it: the list is known no further, so no later span can add to it.
    if (a_spans.m_horizon <= t1) {
      return;
    }
  }
}
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline HitT<T>
levelHit(const LevelView<S>&  a_view,
         const RayT<T>&       a_ray,
         const T              a_tMin,
         const T              a_tMax,
         const RayOptions<T>& a_options) noexcept
{
  return RayForm<LevelView<S>>::firstHit(a_view, a_ray, a_tMin, a_tMax, a_options);
}
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline HitT<T>
levelLeafHit(const LevelView<S>& a_view, const RayT<T>& a_ray, const T a_t, const RayOptions<T>& a_options) noexcept
{
  const auto box   = a_view.computeBoundingVolume();
  const T    width = a_options.m_tolerance * RayDetail::scaleOf(box.getLowCorner(), box.getHighCorner(), a_options);
 
  HitT<T> nearest;
 
  T from = a_t - width;
 
  // As leafHit(): at most three crossings in the window.
  for (int k = 0; k < 3; k++) {
    const HitT<T> hit = RayDetail::levelHit(a_view, a_ray, from, a_t + width, a_options);
 
    if (!hit.hit()) {
      break;
    }
 
    if (!nearest.hit() || std::abs(hit.m_t - a_t) < std::abs(nearest.m_t - a_t)) {
      nearest = hit;
    }
 
    if (hit.m_t >= a_t) {
      break;
    }
 
    from = hit.m_t + T(1e-3) * width;
  }
 
  if (!nearest.hit()) {
    nearest.m_normal = RayDetail::levelNormal(a_view, a_ray.at(a_t), width);
  }
 
  return nearest;
}
 
} // namespace RayDetail
 
template <class T, bool IsAnnulus, class C>
struct RayForm<Fused::Node<Fused::ShiftOp<T, IsAnnulus>, C>, std::enable_if_t<RayDetail::HasLevelForm<C>::value>>
{
  using Node = Fused::Node<Fused::ShiftOp<T, IsAnnulus>, C>;
 
  [[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
  {
    return RayDetail::levelHit(view(a_shape), a_ray, a_tMin, a_tMax, a_options);
  }
 
  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
  {
    RayDetail::levelSpans(view(a_shape), a_ray, a_tMin, a_tMax, a_options, a_firstOnly, a_spans);
  }
 
private:
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline RayDetail::LevelView<C>
  view(const Node& a_shape) noexcept
  {
    return RayDetail::LevelView<C>{RayDetail::fusedOperand<0>(a_shape), a_shape.layout().m_op.m_amount, IsAnnulus};
  }
};
 
} // namespace EBGeometry
 
#endif