EBGeometry_MeshRayForms.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_MESHRAYFORMS_HPP
#define EBGEOMETRY_MESHRAYFORMS_HPP
 
// Std includes
#include <cstddef>
#include <cstdint>
 
// Our includes
#include "EBGeometry_MeshDistanceFunctions.hpp"
#include "EBGeometry_Ray.hpp"
 
namespace EBGeometry {
 
template <class T, size_t K, size_t W>
struct RayForm<TriMeshSDF<T, K, W>>
{
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(const TriMeshSDF<T, K, W>& a_shape,
           const RayT<T>&             a_ray,
           T                          a_tMin,
           T                          a_tMax,
           const RayOptions<T>&       a_options) noexcept;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitStatus
  intersects(const TriMeshSDF<T, K, W>& a_shape,
             const RayT<T>&             a_ray,
             T                          a_tMin,
             T                          a_tMax,
             const RayOptions<T>&       a_options) noexcept;
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const TriMeshSDF<T, K, W>& 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;
};
 
template <class T, size_t K>
struct RayForm<MeshSDF<T, K>>
{
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(
    const MeshSDF<T, K>& a_shape, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitStatus
  intersects(
    const MeshSDF<T, K>& a_shape, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept;
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const MeshSDF<T, 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;
};
 
template <class T>
struct RayForm<FlatMeshSDF<T>>
{
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  static inline HitT<T>
  firstHit(
    const FlatMeshSDF<T>& a_shape, const RayT<T>& a_ray, T a_tMin, T a_tMax, const RayOptions<T>& a_options) noexcept;
 
  EBGEOMETRY_HOST_DEVICE
  static inline void
  spans(const FlatMeshSDF<T>& 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;
};
 
namespace RayDetail {
 
template <class T>
[[nodiscard]] EBGEOMETRY_HOST_DEVICE
inline HitT<T>
faceHit(const RayT<T>& a_ray, const T a_t, const uint32_t a_face, const Vec3T<T>& a_normal) noexcept
{
  HitT<T> hit;
 
  if (a_face != UINT32_MAX && a_t < Math::Limits<T>::infinity()) {
    hit.m_status    = HitStatus::Hit;
    hit.m_t         = a_t;
    hit.m_point     = a_ray.at(a_t);
    hit.m_primitive = a_face;
    hit.m_normal    = a_normal;
    hit.m_entering  = dot(a_ray.m_direction, a_normal) < T(0);
    hit.m_exact     = true;
  }
 
  return hit;
}
 
template <class T>
struct NearestCrossings
{
  static constexpr uint32_t Capacity = 2U * RaySpansT<T>::Capacity + 1U;
 
  struct Crossing
  {
    T m_t;
 
    bool m_entering;
  };
 
  Crossing m_crossings[Capacity];
 
  uint32_t m_count = 0;
 
  T m_tMax;
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  insert(const T a_t, const bool a_entering) noexcept
  {
    if (m_count == Capacity) {
      if (!(a_t < m_crossings[Capacity - 1U].m_t)) {
        return;
      }
 
      m_count--;
    }
 
    uint32_t i = m_count;
 
    while (i > 0 && m_crossings[i - 1U].m_t > a_t) {
      m_crossings[i] = m_crossings[i - 1U];
 
      i--;
    }
 
    m_crossings[i] = Crossing{a_t, a_entering};
 
    m_count++;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  inline T
  horizon() const noexcept
  {
    return (m_count == Capacity) ? m_crossings[Capacity - 1U].m_t : m_tMax;
  }
 
  EBGEOMETRY_HOST_DEVICE
  inline void
  toSpans(const T a_t0, const bool a_inside, RaySpansT<T>& a_spans) const noexcept
  {
    const T horizon = this->horizon();
 
    a_spans.reset(horizon);
 
    if (m_count == 0U) {
      if (a_inside) {
        a_spans.add(RayDetail::openAt(a_t0), RayDetail::openAt(horizon));
      }
 
      return;
    }
 
    CrossingPairer<T> pairer;
 
    for (uint32_t i = 0; i < m_count; i++) {
      const Crossing& c = m_crossings[i];
 
      pairer.take(SpanEndT<T>{c.m_t, 0U, SpanEndT<T>::Exact}, c.m_entering, a_t0, a_spans);
    }
 
    pairer.finish(horizon, a_spans);
  }
};
 
} // namespace RayDetail
 
template <class T, size_t K, size_t W>
EBGEOMETRY_HOST_DEVICE
inline HitT<T>
RayForm<TriMeshSDF<T, K, W>>::firstHit(const TriMeshSDF<T, K, W>& a_shape,
                                       const RayT<T>&             a_ray,
                                       const T                    a_tMin,
                                       const T                    a_tMax,
                                       const RayOptions<T>& /*a_options*/) noexcept
{
  // The state is the nearest crossing so far; it starts at the far end, so that the walk's bound (the
  // state's distance) leaves out every box beyond it.
  struct Nearest
  {
    T        m_t;
    uint32_t m_face;
    Vec3T<T> m_normal;
  };
 
  Nearest nearest{a_tMax, UINT32_MAX, Vec3T<T>::zeros()};
 
  const auto bvh    = a_shape.getRoot();
  const auto groups = bvh.getPrimitives();
 
  const auto leaf = [&groups, &a_ray, a_tMin](Nearest& a_state, size_t a_offset, size_t a_count) noexcept {
    for (size_t i = a_offset; i < a_offset + a_count; i++) {
      uint32_t face   = UINT32_MAX;
      Vec3T<T> normal = Vec3T<T>::zeros();
 
      const T t = groups[static_cast<uint32_t>(i)].rayCrossing(
        a_ray.m_origin, a_ray.m_direction, a_tMin, a_state.m_t, face, normal);
 
      if (face != UINT32_MAX && t <= a_state.m_t) {
        a_state = Nearest{t, face, normal};
      }
    }
  };
 
  const auto bound = [](const Nearest& a_state) noexcept -> T { return a_state.m_t; };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), nearest, leaf, bound);
 
  HitT<T> hit;
 
  if (nearest.m_face != UINT32_MAX) {
    hit.m_status    = HitStatus::Hit;
    hit.m_t         = nearest.m_t;
    hit.m_point     = a_ray.at(nearest.m_t);
    hit.m_primitive = nearest.m_face;
    hit.m_normal    = nearest.m_normal;
    hit.m_entering  = dot(a_ray.m_direction, nearest.m_normal) < T(0);
    hit.m_exact     = true;
  }
 
  return hit;
}
 
template <class T, size_t K, size_t W>
EBGEOMETRY_HOST_DEVICE
inline HitStatus
RayForm<TriMeshSDF<T, K, W>>::intersects(const TriMeshSDF<T, K, W>& a_shape,
                                         const RayT<T>&             a_ray,
                                         const T                    a_tMin,
                                         const T                    a_tMax,
                                         const RayOptions<T>& /*a_options*/) noexcept
{
  // The walk ends once a crossing is found: the bound drops to -infinity, which no box is within.
  bool found = false;
 
  const auto bvh    = a_shape.getRoot();
  const auto groups = bvh.getPrimitives();
 
  const auto leaf = [&groups, &a_ray, a_tMin, a_tMax](bool& a_found, size_t a_offset, size_t a_count) noexcept {
    for (size_t i = a_offset; i < a_offset + a_count && !a_found; i++) {
      uint32_t face   = UINT32_MAX;
      Vec3T<T> normal = Vec3T<T>::zeros();
 
      const T t =
        groups[static_cast<uint32_t>(i)].rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, a_tMax, face, normal);
 
      a_found = face != UINT32_MAX && t <= a_tMax;
    }
  };
 
  const auto bound = [a_tMax](const bool& a_found) noexcept -> T {
    return a_found ? -Math::Limits<T>::infinity() : a_tMax;
  };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), found, leaf, bound);
 
  return found ? HitStatus::Hit : HitStatus::Miss;
}
 
template <class T, size_t K>
EBGEOMETRY_HOST_DEVICE
inline HitT<T>
RayForm<MeshSDF<T, K>>::firstHit(const MeshSDF<T, K>& a_shape,
                                 const RayT<T>&       a_ray,
                                 const T              a_tMin,
                                 const T              a_tMax,
                                 const RayOptions<T>& /*a_options*/) noexcept
{
  // The state is the nearest crossing so far and its face; it starts at the far end, so that the walk's
  // bound leaves out every box beyond it.
  struct Nearest
  {
    T        m_t;
    uint32_t m_face;
  };
 
  Nearest nearest{a_tMax, UINT32_MAX};
 
  const auto bvh     = a_shape.getRoot();
  const auto mesh    = a_shape.getMesh();
  const auto faceIds = bvh.getPrimitives();
 
  const auto leaf = [&faceIds, &mesh, &a_ray, a_tMin](Nearest& a_state, size_t a_offset, size_t a_count) noexcept {
    for (size_t i = a_offset; i < a_offset + a_count; i++) {
      const uint32_t face = faceIds[static_cast<uint32_t>(i)];
      const T        t = mesh.getFace(face).rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, a_state.m_t, mesh);
 
      // A crossing (finite: none is +infinity, however far the interval goes) is within the nearest so far.
      if (t < Math::Limits<T>::infinity()) {
        a_state = Nearest{t, face};
      }
    }
  };
 
  const auto bound = [](const Nearest& a_state) noexcept -> T { return a_state.m_t; };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), nearest, leaf, bound);
 
  return nearest.m_face == UINT32_MAX
           ? HitT<T>()
           : RayDetail::faceHit(a_ray, nearest.m_t, nearest.m_face, mesh.getFace(nearest.m_face).getNormal());
}
 
template <class T, size_t K>
EBGEOMETRY_HOST_DEVICE
inline HitStatus
RayForm<MeshSDF<T, K>>::intersects(const MeshSDF<T, K>& a_shape,
                                   const RayT<T>&       a_ray,
                                   const T              a_tMin,
                                   const T              a_tMax,
                                   const RayOptions<T>& /*a_options*/) noexcept
{
  // The walk ends once a crossing is found: the bound drops to -infinity, which no box is within.
  bool found = false;
 
  const auto bvh     = a_shape.getRoot();
  const auto mesh    = a_shape.getMesh();
  const auto faceIds = bvh.getPrimitives();
 
  const auto leaf = [&faceIds, &mesh, &a_ray, a_tMin, a_tMax](bool& a_found, size_t a_offset, size_t a_count) noexcept {
    for (size_t i = a_offset; i < a_offset + a_count && !a_found; i++) {
      const auto& face = mesh.getFace(faceIds[static_cast<uint32_t>(i)]);
 
      // None is +infinity, which an unbounded interval would contain.
      a_found = face.rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, a_tMax, mesh) < Math::Limits<T>::infinity();
    }
  };
 
  const auto bound = [a_tMax](const bool& a_found) noexcept -> T {
    return a_found ? -Math::Limits<T>::infinity() : a_tMax;
  };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), found, leaf, bound);
 
  return found ? HitStatus::Hit : HitStatus::Miss;
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline HitT<T>
RayForm<FlatMeshSDF<T>>::firstHit(const FlatMeshSDF<T>& a_shape,
                                  const RayT<T>&        a_ray,
                                  const T               a_tMin,
                                  const T               a_tMax,
                                  const RayOptions<T>& /*a_options*/) noexcept
{
  const auto mesh = a_shape.getMesh();
 
  T        nearest = a_tMax;
  uint32_t face    = UINT32_MAX;
 
  for (uint32_t i = 0; i < mesh.numFaces(); i++) {
    const T t = mesh.getFace(i).rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, nearest, mesh);
 
    if (t < Math::Limits<T>::infinity()) {
      nearest = t;
      face    = i;
    }
  }
 
  return face == UINT32_MAX ? HitT<T>() : RayDetail::faceHit(a_ray, nearest, face, mesh.getFace(face).getNormal());
}
 
template <class T, size_t K, size_t W>
EBGEOMETRY_HOST_DEVICE
inline void
RayForm<TriMeshSDF<T, K, W>>::spans(const TriMeshSDF<T, K, W>& a_shape,
                                    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
{
  if (a_firstOnly) {
    RayDetail::firstSpan(a_shape, a_ray, a_tMin, a_tMax, a_options, a_spans);
 
    return;
  }
 
  RayDetail::NearestCrossings<T> crossings;
 
  crossings.m_tMax = a_tMax;
 
  const auto bvh    = a_shape.getRoot();
  const auto groups = bvh.getPrimitives();
 
  const auto leaf =
    [&groups, &a_ray, a_tMin](RayDetail::NearestCrossings<T>& a_state, 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)];
 
        T t[W];
 
        group.rayCrossings(a_ray.m_origin, a_ray.m_direction, a_tMin, a_state.horizon(), t);
 
        for (size_t lane = 0; lane < W; lane++) {
          if (t[lane] < Math::Limits<T>::infinity()) {
            a_state.insert(t[lane], dot(a_ray.m_direction, group.getNormal(lane)) < T(0));
          }
        }
      }
    };
 
  const auto bound = [](const RayDetail::NearestCrossings<T>& a_state) noexcept -> T { return a_state.horizon(); };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), crossings, leaf, bound);
 
  crossings.toSpans(a_tMin, crossings.m_count == 0U && RayDetail::insideAt(a_shape, a_ray.at(a_tMin)), a_spans);
}
 
template <class T, size_t K>
EBGEOMETRY_HOST_DEVICE
inline void
RayForm<MeshSDF<T, K>>::spans(const MeshSDF<T, K>& a_shape,
                              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
{
  if (a_firstOnly) {
    RayDetail::firstSpan(a_shape, a_ray, a_tMin, a_tMax, a_options, a_spans);
 
    return;
  }
 
  RayDetail::NearestCrossings<T> crossings;
 
  crossings.m_tMax = a_tMax;
 
  const auto bvh     = a_shape.getRoot();
  const auto mesh    = a_shape.getMesh();
  const auto faceIds = bvh.getPrimitives();
 
  const auto leaf = [&faceIds, &mesh, &a_ray, a_tMin](
                      RayDetail::NearestCrossings<T>& a_state, size_t a_offset, size_t a_count) noexcept {
    for (size_t i = a_offset; i < a_offset + a_count; i++) {
      const uint32_t face = faceIds[static_cast<uint32_t>(i)];
      const auto&    f    = mesh.getFace(face);
      const T        t    = f.rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, a_state.horizon(), mesh);
 
      if (t < Math::Limits<T>::infinity()) {
        a_state.insert(t, dot(a_ray.m_direction, f.getNormal()) < T(0));
      }
    }
  };
 
  const auto bound = [](const RayDetail::NearestCrossings<T>& a_state) noexcept -> T { return a_state.horizon(); };
 
  bvh.traverse(BVH::RayQuery<T>(a_ray.m_origin, a_ray.m_direction, a_tMin), crossings, leaf, bound);
 
  crossings.toSpans(a_tMin, crossings.m_count == 0U && RayDetail::insideAt(a_shape, a_ray.at(a_tMin)), a_spans);
}
 
template <class T>
EBGEOMETRY_HOST_DEVICE
inline void
RayForm<FlatMeshSDF<T>>::spans(const FlatMeshSDF<T>& a_shape,
                               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
{
  if (a_firstOnly) {
    RayDetail::firstSpan(a_shape, a_ray, a_tMin, a_tMax, a_options, a_spans);
 
    return;
  }
 
  RayDetail::NearestCrossings<T> crossings;
 
  crossings.m_tMax = a_tMax;
 
  const auto mesh = a_shape.getMesh();
 
  for (uint32_t i = 0; i < mesh.numFaces(); i++) {
    const auto& f = mesh.getFace(i);
    const T     t = f.rayCrossing(a_ray.m_origin, a_ray.m_direction, a_tMin, crossings.horizon(), mesh);
 
    if (t < Math::Limits<T>::infinity()) {
      crossings.insert(t, dot(a_ray.m_direction, f.getNormal()) < T(0));
    }
  }
 
  crossings.toSpans(a_tMin, crossings.m_count == 0U && RayDetail::insideAt(a_shape, a_ray.at(a_tMin)), a_spans);
}
 
} // namespace EBGeometry
 
#endif