EBGeometry_RaySpans.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_RAYSPANS_HPP
#define EBGEOMETRY_RAYSPANS_HPP
 
// Std includes
#include <cmath>
#include <cstdint>
 
// Our includes
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_Ray.hpp"
 
namespace EBGeometry {
 
// A span list is left uninitialized past m_count (zeroing it costs local-memory stores on a device), and GCC
// cannot see that only the first m_count spans are ever read.
#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
#endif
 
namespace RayDetail {
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline SpanEndT<T>
openAt(const T a_t) noexcept
{
  return SpanEndT<T>{a_t, 0U, SpanEndT<T>::Open};
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline SpanEndT<T>
endOf(const HitT<T>& a_hit) noexcept
{
  return SpanEndT<T>{a_hit.m_t, 0U, a_hit.m_exact ? SpanEndT<T>::Exact : SpanEndT<T>::Marched};
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
copySpans(const RaySpansT<T>& a_from, RaySpansT<T>& a_to) noexcept
{
  a_to.m_count   = a_from.m_count;
  a_to.m_horizon = a_from.m_horizon;
 
  for (uint32_t i = 0; i < a_from.m_count; i++) {
    a_to.m_spans[i] = a_from.m_spans[i];
  }
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline typename RaySpansT<T>::Span
cut(const typename RaySpansT<T>::Span& a_span, const T a_horizon) noexcept
{
  typename RaySpansT<T>::Span span = a_span;
 
  if (span.m_out.m_t > a_horizon) {
    span.m_out = openAt(a_horizon);
  }
 
  return span;
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
unite(const RaySpansT<T>& a_a, const RaySpansT<T>& a_b, RaySpansT<T>& a_out) noexcept
{
  const T horizon = Math::min(a_a.m_horizon, a_b.m_horizon);
 
  a_out.reset(horizon);
 
  uint32_t i = 0;
  uint32_t j = 0;
 
  bool                        open = false;
  typename RaySpansT<T>::Span current{};
 
  while (true) {
    const bool fromA = i < a_a.m_count && a_a.m_spans[i].m_in.m_t <= horizon;
    const bool fromB = j < a_b.m_count && a_b.m_spans[j].m_in.m_t <= horizon;
 
    if (!fromA && !fromB) {
      break;
    }
 
    // The span that begins first, and of two that begin together the one already inside.
    bool takeA = fromA;
 
    if (fromA && fromB) {
      const SpanEndT<T>& a = a_a.m_spans[i].m_in;
      const SpanEndT<T>& b = a_b.m_spans[j].m_in;
 
      takeA = a.m_t < b.m_t || (a.m_t == b.m_t && !a.crossing());
    }
 
    const typename RaySpansT<T>::Span span = cut(takeA ? a_a.m_spans[i++] : a_b.m_spans[j++], horizon);
 
    if (!open) {
      current = span;
      open    = true;
    }
    else if (span.m_in.m_t <= current.m_out.m_t) {
      if (span.m_out.m_t > current.m_out.m_t || (span.m_out.m_t == current.m_out.m_t && !span.m_out.crossing())) {
        current.m_out = span.m_out;
      }
    }
    else {
      a_out.add(current.m_in, current.m_out);
 
      current = span;
    }
  }
 
  if (open) {
    a_out.add(current.m_in, current.m_out);
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
intersect(const RaySpansT<T>& a_a, const RaySpansT<T>& a_b, RaySpansT<T>& a_out) noexcept
{
  const T horizon = Math::min(a_a.m_horizon, a_b.m_horizon);
 
  a_out.reset(horizon);
 
  uint32_t i = 0;
  uint32_t j = 0;
 
  while (i < a_a.m_count && j < a_b.m_count) {
    if (a_a.m_spans[i].m_in.m_t > horizon || a_b.m_spans[j].m_in.m_t > horizon) {
      break;
    }
 
    const typename RaySpansT<T>::Span a = cut(a_a.m_spans[i], horizon);
    const typename RaySpansT<T>::Span b = cut(a_b.m_spans[j], horizon);
 
    const bool laterA = a.m_in.m_t > b.m_in.m_t || (a.m_in.m_t == b.m_in.m_t && a.m_in.crossing());
    const bool earlyA = a.m_out.m_t < b.m_out.m_t || (a.m_out.m_t == b.m_out.m_t && a.m_out.crossing());
 
    const SpanEndT<T>& in  = laterA ? a.m_in : b.m_in;
    const SpanEndT<T>& out = earlyA ? a.m_out : b.m_out;
 
    if (in.m_t <= out.m_t) {
      a_out.add(in, out);
    }
 
    // The span that ends first can meet nothing more.
    if (a.m_out.m_t <= b.m_out.m_t) {
      i++;
    }
    else {
      j++;
    }
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
cutTo(const T a_horizon, RaySpansT<T>& a_spans) noexcept
{
  uint32_t count = 0;
 
  while (count < a_spans.m_count && a_spans.m_spans[count].m_in.m_t <= a_horizon) {
    if (a_spans.m_spans[count].m_out.m_t > a_horizon) {
      a_spans.m_spans[count].m_out = openAt(a_horizon);
    }
 
    count++;
  }
 
  a_spans.m_count   = count;
  a_spans.m_horizon = a_horizon;
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
keepFirst(RaySpansT<T>& a_spans) noexcept
{
  if (a_spans.m_count == 0U) {
    return;
  }
 
  const typename RaySpansT<T>::Span& first = a_spans.m_spans[0];
 
  if (first.m_in.crossing()) {
    RayDetail::cutTo(first.m_in.m_t, a_spans);
  }
  else if (first.m_out.crossing()) {
    RayDetail::cutTo(first.m_out.m_t, a_spans);
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
uniteInto(RaySpansT<T>& a_spans, const RaySpansT<T>& a_operand) noexcept
{
  const T horizon = Math::min(a_spans.m_horizon, a_operand.m_horizon);
 
  if (horizon < a_spans.m_horizon) {
    RayDetail::cutTo(horizon, a_spans);
  }
 
  // Each of the operand's spans in turn: merged with the list's spans it meets, or inserted between them.
  for (uint32_t k = 0; k < a_operand.m_count && a_operand.m_spans[k].m_in.m_t <= horizon; k++) {
    typename RaySpansT<T>::Span span = RayDetail::cut(a_operand.m_spans[k], horizon);
 
    uint32_t first = 0;
 
    while (first < a_spans.m_count && a_spans.m_spans[first].m_out.m_t < span.m_in.m_t) {
      first++;
    }
 
    // Where two ends meet, the one that is no crossing wins: the ray is inside the other span there.
    uint32_t last = first;
 
    while (last < a_spans.m_count && a_spans.m_spans[last].m_in.m_t <= span.m_out.m_t) {
      const typename RaySpansT<T>::Span& other = a_spans.m_spans[last];
 
      if (other.m_in.m_t < span.m_in.m_t || (other.m_in.m_t == span.m_in.m_t && !other.m_in.crossing())) {
        span.m_in = other.m_in;
      }
      if (other.m_out.m_t > span.m_out.m_t || (other.m_out.m_t == span.m_out.m_t && !other.m_out.crossing())) {
        span.m_out = other.m_out;
      }
 
      last++;
    }
 
    const uint32_t merged = last - first;
 
    if (merged == 0U) {
      // A new span: a full list keeps its first spans, and is known only up to where the dropped one begins.
      if (a_spans.m_count == RaySpansT<T>::Capacity) {
        if (first == a_spans.m_count) {
          a_spans.m_horizon = Math::min(a_spans.m_horizon, span.m_in.m_t);
 
          break;
        }
 
        a_spans.m_horizon = Math::min(a_spans.m_horizon, a_spans.m_spans[a_spans.m_count - 1U].m_in.m_t);
        a_spans.m_count--;
      }
 
      for (uint32_t m = a_spans.m_count; m > first; m--) {
        a_spans.m_spans[m] = a_spans.m_spans[m - 1U];
      }
 
      a_spans.m_count++;
    }
    else if (merged > 1U) {
      for (uint32_t m = last; m < a_spans.m_count; m++) {
        a_spans.m_spans[m - merged + 1U] = a_spans.m_spans[m];
      }
 
      a_spans.m_count -= merged - 1U;
    }
 
    a_spans.m_spans[first] = span;
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
intersectInto(RaySpansT<T>& a_spans, const RaySpansT<T>& a_operand, RaySpansT<T>& a_scratch) noexcept
{
  if (a_operand.m_count == 0U || a_spans.m_count == 0U) {
    a_spans.reset(Math::min(a_spans.m_horizon, a_operand.m_horizon));
 
    return;
  }
 
  RayDetail::intersect(a_spans, a_operand, a_scratch);
  RayDetail::copySpans(a_scratch, a_spans);
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
complement(const RaySpansT<T>& a_a, const T a_t0, RaySpansT<T>& a_out) noexcept
{
  const T horizon = a_a.m_horizon;
 
  a_out.reset(horizon);
 
  // The gap the ray is in: from where it starts, while there is one.
  SpanEndT<T> start = openAt(a_t0);
  bool        gap   = true;
 
  for (uint32_t i = 0; i < a_a.m_count && a_a.m_spans[i].m_in.m_t <= horizon; i++) {
    const typename RaySpansT<T>::Span span = cut(a_a.m_spans[i], horizon);
 
    // A span that begins at the start of the interval, inside, has no gap before it.
    if (span.m_in.crossing() && gap) {
      a_out.add(start, span.m_in);
    }
 
    // A span that ends where nothing more is known leaves no gap after it.
    gap = span.m_out.crossing();
 
    if (!gap) {
      break;
    }
 
    start = span.m_out;
  }
 
  if (gap) {
    a_out.add(start, openAt(horizon));
  }
}
 
template <class T>
struct MappedRay
{
  RayT<T> m_ray;
 
  T m_scale;
 
  T m_t0;
 
  T m_t1;
};
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline MappedRay<T>
mapRay(const T* a_map, const RayT<T>& a_ray, const T a_t0, const T a_t1) noexcept
{
  const T*        m = a_map;
  const Vec3T<T>& o = a_ray.m_origin;
  const Vec3T<T>& d = a_ray.m_direction;
 
  const Vec3T<T> origin(m[0] * o[0] + m[1] * o[1] + m[2] * o[2] + m[9],
                        m[3] * o[0] + m[4] * o[1] + m[5] * o[2] + m[10],
                        m[6] * o[0] + m[7] * o[1] + m[8] * o[2] + m[11]);
  const Vec3T<T> direction(m[0] * d[0] + m[1] * d[1] + m[2] * d[2],
                           m[3] * d[0] + m[4] * d[1] + m[5] * d[2],
                           m[6] * d[0] + m[7] * d[1] + m[8] * d[2]);
 
  const T scale = direction.length();
 
  return MappedRay<T>{RayT<T>{origin, direction / scale}, scale, a_t0 * scale, a_t1 * scale};
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
mapBack(const MappedRay<T>& a_mapped, const T a_t0, const T a_t1, RaySpansT<T>& a_spans) noexcept
{
  const auto distance = [&a_mapped, a_t0, a_t1](const T a_t) noexcept -> T {
    if (a_t <= a_mapped.m_t0) {
      return a_t0;
    }
    if (a_t >= a_mapped.m_t1) {
      return a_t1;
    }
 
    // Rounded, the quotient can land just outside the original interval.
    return Math::clamp(a_t / a_mapped.m_scale, a_t0, a_t1);
  };
 
  a_spans.m_horizon = distance(a_spans.m_horizon);
 
  for (uint32_t i = 0; i < a_spans.m_count; i++) {
    a_spans.m_spans[i].m_in.m_t  = distance(a_spans.m_spans[i].m_in.m_t);
    a_spans.m_spans[i].m_out.m_t = distance(a_spans.m_spans[i].m_out.m_t);
  }
}
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline Vec3T<T>
mapNormalBack(const T* a_map, const Vec3T<T>& a_normal) noexcept
{
  const T*        m = a_map;
  const Vec3T<T>& n = a_normal;
  const Vec3T<T>  g(m[0] * n[0] + m[3] * n[1] + m[6] * n[2],
                   m[1] * n[0] + m[4] * n[1] + m[7] * n[2],
                   m[2] * n[0] + m[5] * n[1] + m[8] * n[2]);
  const T         length = g.length();
 
  return (length > T(0)) ? g / length : g;
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
nameLeaf(const uint32_t a_leaf, RaySpansT<T>& a_spans) noexcept
{
  for (uint32_t i = 0; i < a_spans.m_count; i++) {
    a_spans.m_spans[i].m_in.m_leaf  = a_leaf;
    a_spans.m_spans[i].m_out.m_leaf = a_leaf;
  }
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
offsetLeaves(const uint32_t a_offset, RaySpansT<T>& a_spans) noexcept
{
  for (uint32_t i = 0; i < a_spans.m_count; i++) {
    a_spans.m_spans[i].m_in.m_leaf += a_offset;
    a_spans.m_spans[i].m_out.m_leaf += a_offset;
  }
}
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline bool
insideAt(const S& a_shape, const Vec3T<T>& a_point) noexcept
{
  const auto      box = a_shape.computeBoundingVolume();
  const Vec3T<T>& lo  = box.getLowCorner();
  const Vec3T<T>& hi  = box.getHighCorner();
 
  for (size_t dir = 0; dir < 3; dir++) {
    if (!(a_point[dir] >= lo[dir] && a_point[dir] <= hi[dir])) {
      return false;
    }
  }
 
  return RayDetail::valueAt(a_shape, a_point) < T(0);
}
 
template <class T, class S>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
EBGEOMETRY_DEVICE_NOINLINE
inline HitT<T>
shapeHit(const S& a_shape, const RayT<T>& a_ray, const T a_t0, const T a_t1, const RayOptions<T>& a_options) noexcept
{
  return EBGeometry::firstHit(a_shape, a_ray, a_t0, a_t1, a_options);
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
inline void
firstSpan(const S&             a_shape,
          const RayT<T>&       a_ray,
          const T              a_t0,
          const T              a_t1,
          const RayOptions<T>& a_options,
          RaySpansT<T>&        a_spans) noexcept
{
  a_spans.reset(a_t1);
 
  const HitT<T> hit = RayDetail::shapeHit(a_shape, a_ray, a_t0, a_t1, a_options);
 
  switch (hit.m_status) {
  case HitStatus::Hit: {
    if (hit.m_entering) {
      a_spans.add(endOf(hit), openAt(hit.m_t));
    }
    else {
      a_spans.add(openAt(a_t0), endOf(hit));
    }
 
    a_spans.m_horizon = hit.m_t;
 
    break;
  }
  case HitStatus::Unresolved: {
    a_spans.m_horizon = Math::max(a_t0, Math::min(hit.m_t, a_t1));
 
    break;
  }
  default: {
    if (RayDetail::insideAt(a_shape, a_ray.at(a_t0))) {
      a_spans.add(openAt(a_t0), openAt(a_t1));
    }
 
    break;
  }
  }
}
 
template <class T>
struct CrossingPairer
{
  SpanEndT<T> m_start;
 
  bool m_inside = false;
 
  bool m_first = true;
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  take(const SpanEndT<T>& a_end, const bool a_entering, const T a_t0, RaySpansT<T>& a_spans) noexcept
  {
    if (a_entering) {
      if (!m_inside) {
        m_start  = a_end;
        m_inside = true;
      }
    }
    else if (m_inside) {
      a_spans.add(m_start, a_end);
 
      m_inside = false;
    }
    else if (m_first) {
      a_spans.add(openAt(a_t0), a_end);
    }
 
    m_first = false;
  }
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  finish(const T a_t, RaySpansT<T>& a_spans) noexcept
  {
    if (m_inside) {
      a_spans.add(m_start, openAt(a_t));
 
      m_inside = false;
    }
  }
};
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
inline void
crossingSpans(const S&             a_shape,
              const RayT<T>&       a_ray,
              const T              a_t0,
              const T              a_t1,
              const RayOptions<T>& a_options,
              RaySpansT<T>&        a_spans) noexcept
{
  a_spans.reset(a_t1);
 
  const auto box   = a_shape.computeBoundingVolume();
  const T    scale = RayDetail::scaleOf(box.getLowCorner(), box.getHighCorner(), a_options);
  const T    step  = a_options.m_tolerance * scale;
 
  CrossingPairer<T> pairer;
 
  T t = a_t0;
 
  // Two crossings per span, and one more to find out the list is full.
  for (uint32_t k = 0; k < 2U * RaySpansT<T>::Capacity + 1U; k++) {
    const HitT<T> hit = RayDetail::shapeHit(a_shape, a_ray, t, a_t1, a_options);
 
    if (hit.m_status == HitStatus::Miss) {
      break;
    }
 
    if (hit.m_status == HitStatus::Unresolved) {
      const T horizon = Math::max(a_t0, Math::min(hit.m_t, a_t1));
 
      pairer.finish(horizon, a_spans);
 
      a_spans.m_horizon = Math::min(a_spans.m_horizon, horizon);
 
      return;
    }
 
    pairer.take(endOf(hit), hit.m_entering, a_t0, a_spans);
 
    if (a_spans.m_horizon < a_t1) {
      return;
    }
 
    t = hit.m_t + step;
 
    if (!(t <= a_t1)) {
      break;
    }
 
    if (k + 1U == 2U * RaySpansT<T>::Capacity + 1U) {
      // Out of searches with the interval not done: known only up to here.
      pairer.finish(t, a_spans);
 
      a_spans.m_horizon = Math::min(a_spans.m_horizon, t);
 
      return;
    }
  }
 
  if (pairer.m_first) {
    if (RayDetail::insideAt(a_shape, a_ray.at(a_t0))) {
      a_spans.add(openAt(a_t0), openAt(a_t1));
    }
  }
  else {
    pairer.finish(a_t1, a_spans);
  }
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
inline void
shapeSpans(const S&             a_shape,
           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
{
  if constexpr (HasRaySpans<S, T>::value) {
    RayForm<S>::spans(a_shape, a_ray, a_t0, a_t1, a_options, a_firstOnly, a_spans);
  }
  else {
    if (a_firstOnly) {
      RayDetail::firstSpan(a_shape, a_ray, a_t0, a_t1, a_options, a_spans);
    }
    else {
      RayDetail::crossingSpans(a_shape, a_ray, a_t0, a_t1, a_options, a_spans);
    }
  }
}
 
template <class T, class S>
EBGEOMETRY_HOST_DEVICE
inline void
spansWithin(const S&             a_shape,
            const RayT<T>&       a_ray,
            const RaySpansT<T>&  a_within,
            const RayOptions<T>& a_options,
            const bool           a_firstOnly,
            RaySpansT<T>&        a_piece,
            RaySpansT<T>&        a_spans) noexcept
{
  a_spans.reset(a_within.m_horizon);
 
  for (uint32_t i = 0; i < a_within.m_count && a_within.m_spans[i].m_in.m_t <= a_within.m_horizon; i++) {
    const T t0 = a_within.m_spans[i].m_in.m_t;
    const T t1 = Math::min(a_within.m_spans[i].m_out.m_t, a_within.m_horizon);
 
    RayDetail::shapeSpans(a_shape, a_ray, t0, t1, a_options, a_firstOnly, a_piece);
 
    for (uint32_t k = 0; k < a_piece.m_count && a_piece.m_spans[k].m_in.m_t <= a_piece.m_horizon; k++) {
      const typename RaySpansT<T>::Span span = RayDetail::cut(a_piece.m_spans[k], a_piece.m_horizon);
 
      a_spans.add(span.m_in, span.m_out);
    }
 
    if (a_piece.m_horizon < t1 || (a_firstOnly && a_piece.m_count > 0U)) {
      a_spans.m_horizon = Math::min(a_spans.m_horizon, a_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
inline HitT<T>
leafHit(const S& a_shape, const RayT<T>& a_ray, const T a_t, const RayOptions<T>& a_options) noexcept
{
  const auto box   = a_shape.computeBoundingVolume();
  const T    delta = a_options.m_tolerance * RayDetail::scaleOf(box.getLowCorner(), box.getHighCorner(), a_options);
 
  HitT<T> nearest;
 
  T from = a_t - delta;
 
  // At most three crossings in the window: a ray grazing an edge may cross two faces there.
  for (int k = 0; k < 3; k++) {
    const HitT<T> hit = RayDetail::shapeHit(a_shape, a_ray, from, a_t + delta, 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) * delta;
  }
 
  if (!nearest.hit()) {
    const Vec3T<T> point = a_ray.at(a_t);
 
    nearest.m_normal    = RayDetail::gradientNormal(a_shape, point, delta);
    nearest.m_primitive = RayDetail::primitiveOf(a_shape, point);
  }
 
  return nearest;
}
 
template <class T, class Trace, class Resolve>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline HitT<T>
traceSpans(
  const Trace& a_trace, const Resolve& a_resolve, const RayT<T>& a_ray, const T a_tMin, const T a_tMax) noexcept
{
  constexpr uint32_t MaxPasses = 64;
 
  HitT<T> hit;
 
  if (!(a_tMin <= a_tMax)) {
    return hit;
  }
 
  T    t0  = a_tMin;
  bool all = false;
 
  RaySpansT<T> spans;
 
  for (uint32_t pass = 0; pass < MaxPasses; pass++) {
    a_trace(t0, a_tMax, all, spans);
 
    const T horizon = spans.m_horizon;
 
    for (uint32_t i = 0; i < spans.m_count; i++) {
      for (int side = 0; side < 2; side++) {
        const SpanEndT<T>& end = (side == 0) ? spans.m_spans[i].m_in : spans.m_spans[i].m_out;
 
        if (end.crossing() && end.m_t >= t0 && end.m_t <= horizon) {
          const HitT<T> leaf = a_resolve(end.m_t, end.m_leaf);
 
          hit.m_status    = HitStatus::Hit;
          hit.m_t         = end.m_t;
          hit.m_point     = a_ray.at(end.m_t);
          hit.m_primitive = leaf.m_primitive;
          hit.m_normal    = leaf.m_normal;
          hit.m_entering  = (side == 0);
          hit.m_exact     = (end.m_kind == SpanEndT<T>::Exact);
 
          // The leaf's normal faces out of the leaf; the composition's faces the way the ray crosses. For a ray
          // along the surface the dot product is rounding noise, and either orientation may come out.
          if ((dot(hit.m_normal, a_ray.m_direction) < T(0)) != hit.m_entering) {
            hit.m_normal = -hit.m_normal;
          }
 
          return hit;
        }
      }
    }
 
    // No crossing up to the horizon: none at all if that is the end, otherwise start again from it.
    if (!(horizon < a_tMax)) {
      return hit;
    }
 
    if (horizon > t0) {
      t0  = horizon;
      all = false;
    }
    else if (!all) {
      all = true;
    }
    else {
      break;
    }
  }
 
  hit.m_status = HitStatus::Unresolved;
  hit.m_t      = t0;
  hit.m_point  = a_ray.at(t0);
 
  return hit;
}
 
} // namespace RayDetail
 
#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC diagnostic pop
#endif
 
} // namespace EBGeometry
 
#endif