EBGeometry_AnalyticDistanceFunctions.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2022 Robert Marskar <robert.marskar@sintef.no>
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef EBGEOMETRY_ANALYTICDISTANCEFUNCTIONS_HPP
#define EBGEOMETRY_ANALYTICDISTANCEFUNCTIONS_HPP
 
// Std includes
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <random>
#include <type_traits>
 
// Our includes
#include "EBGeometry_Array.hpp"
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_Constants.hpp"
#include "EBGeometry_DistanceQuality.hpp"
#include "EBGeometry_GPU.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_Math.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
template <class T>
class PlaneSDF
{
  static_assert(std::is_floating_point_v<T>, "PlaneSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  PlaneSDF() = default;
 
  PlaneSDF(const Vec3T<T>& a_point, const Vec3T<T>& a_normal) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_point[0]) && std::isfinite(a_point[1]) && std::isfinite(a_point[2]),
                       "PlaneSDF: the point must be finite (%g, %g, %g)",
                       double(a_point[0]),
                       double(a_point[1]),
                       double(a_point[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_normal[0]) && std::isfinite(a_normal[1]) && std::isfinite(a_normal[2]),
                       "PlaneSDF: the normal must be finite (%g, %g, %g)",
                       double(a_normal[0]),
                       double(a_normal[1]),
                       double(a_normal[2]));
    EBGEOMETRY_REQUIRE(a_normal.length() > T(0),
                       "PlaneSDF: the normal must be nonzero (%g, %g, %g)",
                       double(a_normal[0]),
                       double(a_normal[1]),
                       double(a_normal[2]));
 
    m_point  = a_point;
    m_normal = a_normal / a_normal.length();
  }
 
  PlaneSDF(const PlaneSDF&) = default;
 
  PlaneSDF(PlaneSDF&&) = default;
 
  PlaneSDF&
  operator=(const PlaneSDF&) = default;
 
  PlaneSDF&
  operator=(PlaneSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
    EBGEOMETRY_EXPECT(std::abs(m_normal.length() - T(1)) < std::sqrt(Math::Limits<T>::epsilon()));
 
    return dot((a_point - m_point), m_normal);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return BoundingVolumes::AABBT<T>(Vec3T<T>::min(), Vec3T<T>::max());
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getPoint() const noexcept
  {
    return m_point;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getNormal() const noexcept
  {
    return m_normal;
  }
 
protected:
  Vec3T<T> m_point = Vec3T<T>::zeros();
 
  Vec3T<T> m_normal = Vec3T<T>(T(0), T(1), T(0));
};
 
template <class T>
class SphereSDF
{
  static_assert(std::is_floating_point_v<T>, "SphereSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  SphereSDF() = default;
 
  SphereSDF(const Vec3T<T>& a_center, const T& a_radius) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center[0]) && std::isfinite(a_center[1]) && std::isfinite(a_center[2]),
                       "SphereSDF: the center must be finite (%g, %g, %g)",
                       double(a_center[0]),
                       double(a_center[1]),
                       double(a_center[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_radius) && a_radius > T(0),
                       "SphereSDF: the radius must be finite and positive (%g)",
                       double(a_radius));
 
    m_center = a_center;
    m_radius = a_radius;
  }
 
  SphereSDF(const SphereSDF&) = default;
 
  SphereSDF(SphereSDF&&) = default;
 
  SphereSDF&
  operator=(const SphereSDF&) = default;
 
  SphereSDF&
  operator=(SphereSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  EBGEOMETRY_HOST_DEVICE
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  EBGEOMETRY_HOST_DEVICE
  T&
  getRadius() noexcept
  {
    return m_radius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    return (a_point - m_center).length() - m_radius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return BoundingVolumes::AABBT<T>(m_center - m_radius * Vec3T<T>::ones(), m_center + m_radius * Vec3T<T>::ones());
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_radius = T(1);
};
 
template <class T>
class BoxSDF
{
  static_assert(std::is_floating_point_v<T>, "BoxSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  BoxSDF() = default;
 
  BoxSDF(const Vec3T<T>& a_loCorner, const Vec3T<T>& a_hiCorner) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_loCorner[0]) && std::isfinite(a_loCorner[1]) && std::isfinite(a_loCorner[2]),
                       "BoxSDF: the low corner must be finite (%g, %g, %g)",
                       double(a_loCorner[0]),
                       double(a_loCorner[1]),
                       double(a_loCorner[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_hiCorner[0]) && std::isfinite(a_hiCorner[1]) && std::isfinite(a_hiCorner[2]),
                       "BoxSDF: the high corner must be finite (%g, %g, %g)",
                       double(a_hiCorner[0]),
                       double(a_hiCorner[1]),
                       double(a_hiCorner[2]));
    EBGEOMETRY_REQUIRE(
      a_loCorner[0] < a_hiCorner[0] && a_loCorner[1] < a_hiCorner[1] && a_loCorner[2] < a_hiCorner[2],
      "BoxSDF: the low corner (%g, %g, %g) must be below the high corner (%g, %g, %g) in every direction",
      double(a_loCorner[0]),
      double(a_loCorner[1]),
      double(a_loCorner[2]),
      double(a_hiCorner[0]),
      double(a_hiCorner[1]),
      double(a_hiCorner[2]));
 
    m_loCorner = a_loCorner;
    m_hiCorner = a_hiCorner;
  }
 
  BoxSDF(const BoxSDF&) = default;
 
  BoxSDF(BoxSDF&&) = default;
 
  BoxSDF&
  operator=(const BoxSDF&) = default;
 
  BoxSDF&
  operator=(BoxSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getLowCorner() const noexcept
  {
    return m_loCorner;
  }
 
  EBGEOMETRY_HOST_DEVICE
  Vec3T<T>&
  getLowCorner() noexcept
  {
    return m_loCorner;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getHighCorner() const noexcept
  {
    return m_hiCorner;
  }
 
  EBGEOMETRY_HOST_DEVICE
  Vec3T<T>&
  getHighCorner() noexcept
  {
    return m_hiCorner;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    // For each coordinate direction, we have delta[dir] if a_point[dir] falls
    // between xLo and xHi. In this case delta[dir] will be the signed distance
    // to the closest box face in the dir-direction. Otherwise, if a_point[dir]
    // is outside the corner we have delta[dir] > 0.
    const Vec3T<T> delta(Math::max(m_loCorner[0] - a_point[0], a_point[0] - m_hiCorner[0]),
                         Math::max(m_loCorner[1] - a_point[1], a_point[1] - m_hiCorner[1]),
                         Math::max(m_loCorner[2] - a_point[2], a_point[2] - m_hiCorner[2]));
 
    // Note: max is max(Vec3T<T>, Vec3T<T>) and not std::max. It returns a
    // vector with coordinate-wise largest components. Note that the first part
    // Math::min(...) is the signed distance on the inside of the box (delta will
    // have negative components). The other part max(Vec3T<T>::zeros(), ...) is
    // for outside the box.
    const T d = Math::min(T(0.0), delta[delta.maxDir(false)]) + max(Vec3T<T>::zeros(), delta).length();
 
    return d;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return BoundingVolumes::AABBT<T>(m_loCorner, m_hiCorner);
  }
 
protected:
  Vec3T<T> m_loCorner = Vec3T<T>(T(-0.5), T(-0.5), T(-0.5));
 
  Vec3T<T> m_hiCorner = Vec3T<T>(T(0.5), T(0.5), T(0.5));
};
 
template <class T>
class TorusSDF
{
  static_assert(std::is_floating_point_v<T>, "TorusSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  TorusSDF() = default;
 
  TorusSDF(const Vec3T<T>& a_center, const T& a_majorRadius, const T& a_minorRadius) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center[0]) && std::isfinite(a_center[1]) && std::isfinite(a_center[2]),
                       "TorusSDF: the center must be finite (%g, %g, %g)",
                       double(a_center[0]),
                       double(a_center[1]),
                       double(a_center[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_majorRadius) && a_majorRadius > T(0),
                       "TorusSDF: the major radius must be finite and positive (%g)",
                       double(a_majorRadius));
    EBGEOMETRY_REQUIRE(std::isfinite(a_minorRadius) && a_minorRadius > T(0),
                       "TorusSDF: the minor radius must be finite and positive (%g)",
                       double(a_minorRadius));
    EBGEOMETRY_REQUIRE(a_minorRadius < a_majorRadius,
                       "TorusSDF: the minor radius (%g) must be less than the major radius (%g)",
                       double(a_minorRadius),
                       double(a_majorRadius));
 
    m_center      = a_center;
    m_majorRadius = a_majorRadius;
    m_minorRadius = a_minorRadius;
  }
 
  TorusSDF(const TorusSDF&) = default;
 
  TorusSDF(TorusSDF&&) = default;
 
  TorusSDF&
  operator=(const TorusSDF&) = default;
 
  TorusSDF&
  operator=(TorusSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  EBGEOMETRY_HOST_DEVICE
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getMajorRadius() const noexcept
  {
    return m_majorRadius;
  }
 
  EBGEOMETRY_HOST_DEVICE
  T&
  getMajorRadius() noexcept
  {
    return m_majorRadius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getMinorRadius() const noexcept
  {
    return m_minorRadius;
  }
 
  EBGEOMETRY_HOST_DEVICE
  T&
  getMinorRadius() noexcept
  {
    return m_minorRadius;
  }
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    const Vec3T<T> p   = a_point - m_center;
    const T        rho = std::sqrt(p[0] * p[0] + p[2] * p[2]) - m_majorRadius;
    const T        d   = std::sqrt(rho * rho + p[1] * p[1]) - m_minorRadius;
 
    return d;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    const T        outer = m_majorRadius + m_minorRadius;
    const Vec3T<T> half(outer, m_minorRadius, outer);
 
    return BoundingVolumes::AABBT<T>(m_center - half, m_center + half);
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_majorRadius = T(1);
 
  T m_minorRadius = T(0.5);
};
 
template <class T>
class CylinderSDF
{
  static_assert(std::is_floating_point_v<T>, "CylinderSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  CylinderSDF() = default;
 
  CylinderSDF(const Vec3T<T>& a_center1, const Vec3T<T>& a_center2, const T& a_radius) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center1[0]) && std::isfinite(a_center1[1]) && std::isfinite(a_center1[2]),
                       "CylinderSDF: the first center must be finite (%g, %g, %g)",
                       double(a_center1[0]),
                       double(a_center1[1]),
                       double(a_center1[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_center2[0]) && std::isfinite(a_center2[1]) && std::isfinite(a_center2[2]),
                       "CylinderSDF: the second center must be finite (%g, %g, %g)",
                       double(a_center2[0]),
                       double(a_center2[1]),
                       double(a_center2[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_radius) && a_radius > T(0),
                       "CylinderSDF: the radius must be finite and positive (%g)",
                       double(a_radius));
    EBGEOMETRY_REQUIRE((a_center2 - a_center1).length() > T(0),
                       "CylinderSDF: the two centers must differ, but both are (%g, %g, %g)",
                       double(a_center1[0]),
                       double(a_center1[1]),
                       double(a_center1[2]));
 
    m_center1 = a_center1;
    m_center2 = a_center2;
    m_radius  = a_radius;
 
    m_center = (m_center2 + m_center1) * T(0.5);
    m_length = (m_center2 - m_center1).length();
    m_axis   = (m_center2 - m_center1) / m_length;
  }
 
  CylinderSDF(const CylinderSDF&) = default;
 
  CylinderSDF(CylinderSDF&&) = default;
 
  CylinderSDF&
  operator=(const CylinderSDF&) = default;
 
  CylinderSDF&
  operator=(CylinderSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter1() const noexcept
  {
    return m_center1;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter2() const noexcept
  {
    return m_center2;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    T d = Math::Limits<T>::infinity();
 
    if (m_length > T(0.0) && m_radius > T(0.0)) {
      EBGEOMETRY_EXPECT(std::abs(m_axis.length() - T(1)) < std::sqrt(Math::Limits<T>::epsilon()));
 
      const Vec3T<T> point = a_point - m_center;
      const T        para  = dot(point, m_axis);
      const Vec3T<T> ortho = point - para * m_axis; // Distance from cylinder axis.
 
      // w: Distance from cylinder wall. < 0 on inside and > 0 on outside.
      // h: Distance from cylinder top.  < 0 on inside and > 0 on outside.
      const T w = ortho.length() - m_radius;
      const T h = std::abs(para) - T(0.5) * m_length;
 
      constexpr T zero = T(0.0);
 
      if (w <= zero && h <= zero) {
        // Inside cylinder
        d = (std::abs(w) < std::abs(h)) ? w : h;
      }
      else if (w <= zero && h > zero) {
        // Above one of the endcaps.
        d = h;
      }
      else if (w > zero && h < zero) {
        // Outside radius but between the endcaps.
        d = w;
      }
      else {
        d = std::sqrt(w * w + h * h);
      }
    }
 
    return d;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    Vec3T<T> reach;
 
    for (size_t dir = 0; dir < 3; dir++) {
      reach[dir] = m_radius * std::sqrt(Math::max(T(0), T(1) - m_axis[dir] * m_axis[dir]));
    }
 
    return BoundingVolumes::AABBT<T>(min(m_center1, m_center2) - reach, max(m_center1, m_center2) + reach);
  }
 
protected:
  Vec3T<T> m_center1 = Vec3T<T>(T(0), T(-0.5), T(0));
 
  Vec3T<T> m_center2 = Vec3T<T>(T(0), T(0.5), T(0));
 
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  Vec3T<T> m_axis = Vec3T<T>(T(0), T(1), T(0));
 
  T m_length = T(1);
 
  T m_radius = T(1);
};
 
template <class T>
class InfiniteCylinderSDF
{
  static_assert(std::is_floating_point_v<T>, "InfiniteCylinderSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  InfiniteCylinderSDF() = default;
 
  InfiniteCylinderSDF(const Vec3T<T>& a_center, const T& a_radius, const size_t a_axis) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center[0]) && std::isfinite(a_center[1]) && std::isfinite(a_center[2]),
                       "InfiniteCylinderSDF: the center must be finite (%g, %g, %g)",
                       double(a_center[0]),
                       double(a_center[1]),
                       double(a_center[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_radius) && a_radius > T(0),
                       "InfiniteCylinderSDF: the radius must be finite and positive (%g)",
                       double(a_radius));
    EBGEOMETRY_REQUIRE(a_axis < 3U, "InfiniteCylinderSDF: the axis must be 0, 1 or 2 (%zu)", a_axis);
 
    m_center = a_center;
    m_radius = a_radius;
    m_axis   = a_axis;
  }
 
  InfiniteCylinderSDF(const InfiniteCylinderSDF&) = default;
 
  InfiniteCylinderSDF(InfiniteCylinderSDF&&) = default;
 
  InfiniteCylinderSDF&
  operator=(const InfiniteCylinderSDF&) = default;
 
  InfiniteCylinderSDF&
  operator=(InfiniteCylinderSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    Vec3T<T> delta = a_point - m_center;
    delta[m_axis]  = 0.0;
 
    const T d = delta.length() - m_radius;
 
    return d;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    Vec3T<T> lo = m_center - m_radius * Vec3T<T>::ones();
    Vec3T<T> hi = m_center + m_radius * Vec3T<T>::ones();
 
    lo[m_axis] = -Math::Limits<T>::max();
    hi[m_axis] = Math::Limits<T>::max();
 
    return BoundingVolumes::AABBT<T>(lo, hi);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const size_t&
  getAxis() const noexcept
  {
    return m_axis;
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_radius = T(1);
 
  size_t m_axis = 1U;
};
 
template <class T>
class CapsuleSDF
{
  static_assert(std::is_floating_point_v<T>, "CapsuleSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  CapsuleSDF() = default;
 
  CapsuleSDF(const Vec3T<T>& a_tip1, const Vec3T<T>& a_tip2, const T& a_radius) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_tip1[0]) && std::isfinite(a_tip1[1]) && std::isfinite(a_tip1[2]),
                       "CapsuleSDF: the first tip must be finite (%g, %g, %g)",
                       double(a_tip1[0]),
                       double(a_tip1[1]),
                       double(a_tip1[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_tip2[0]) && std::isfinite(a_tip2[1]) && std::isfinite(a_tip2[2]),
                       "CapsuleSDF: the second tip must be finite (%g, %g, %g)",
                       double(a_tip2[0]),
                       double(a_tip2[1]),
                       double(a_tip2[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_radius) && a_radius > T(0),
                       "CapsuleSDF: the radius must be finite and positive (%g)",
                       double(a_radius));
    EBGEOMETRY_REQUIRE((a_tip2 - a_tip1).length() >= T(2) * a_radius,
                       "CapsuleSDF: the tips must be at least two radii apart (distance %g, radius %g)",
                       double((a_tip2 - a_tip1).length()),
                       double(a_radius));
 
    const Vec3T<T> axis = (a_tip2 - a_tip1) / length(a_tip2 - a_tip1);
 
    m_center1 = a_tip1 + a_radius * axis;
    m_center2 = a_tip2 - a_radius * axis;
    m_radius  = a_radius;
  }
 
  CapsuleSDF(const CapsuleSDF&) = default;
 
  CapsuleSDF(CapsuleSDF&&) = default;
 
  CapsuleSDF&
  operator=(const CapsuleSDF&) = default;
 
  CapsuleSDF&
  operator=(CapsuleSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    const Vec3T<T> v1 = a_point - m_center1;
    const Vec3T<T> v2 = m_center2 - m_center1;
 
    // Tips exactly two radii apart put both hemisphere centers at the same point: a sphere. Guard the
    // projection, which would otherwise divide zero by zero.
    const T len2 = dot(v2, v2);
 
    const T h = (len2 > T(0)) ? Math::clamp(dot(v1, v2) / len2, T(0.0), T(1.0)) : T(0.0);
    const T d = length(v1 - h * v2) - m_radius;
 
    return d;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return BoundingVolumes::AABBT<T>(min(m_center1, m_center2) - m_radius * Vec3T<T>::ones(),
                                     max(m_center1, m_center2) + m_radius * Vec3T<T>::ones());
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter1() const noexcept
  {
    return m_center1;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const Vec3T<T>&
  getCenter2() const noexcept
  {
    return m_center2;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
protected:
  Vec3T<T> m_center1 = Vec3T<T>(T(0), T(-0.5), T(0));
 
  Vec3T<T> m_center2 = Vec3T<T>(T(0), T(0.5), T(0));
 
  T m_radius = T(0.5);
};
 
template <class T>
class InfiniteConeSDF
{
  static_assert(std::is_floating_point_v<T>, "InfiniteConeSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  InfiniteConeSDF() = default;
 
  InfiniteConeSDF(const Vec3T<T>& a_tip, const T& a_angle) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_tip[0]) && std::isfinite(a_tip[1]) && std::isfinite(a_tip[2]),
                       "InfiniteConeSDF: the tip must be finite (%g, %g, %g)",
                       double(a_tip[0]),
                       double(a_tip[1]),
                       double(a_tip[2]));
    EBGEOMETRY_REQUIRE(a_angle > T(0) && a_angle < T(180),
                       "InfiniteConeSDF: the angle must be in (0, 180) degrees (%g)",
                       double(a_angle));
 
    m_tip  = a_tip;
    m_c[0] = std::sin(T(0.5) * a_angle * pi<T> / T(180));
    m_c[1] = std::cos(T(0.5) * a_angle * pi<T> / T(180));
  }
 
  InfiniteConeSDF(const InfiniteConeSDF&) = default;
 
  InfiniteConeSDF(InfiniteConeSDF&&) = default;
 
  InfiniteConeSDF&
  operator=(const InfiniteConeSDF&) = default;
 
  InfiniteConeSDF&
  operator=(InfiniteConeSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
    EBGEOMETRY_EXPECT(std::abs(length(m_c) - T(1)) < std::sqrt(Math::Limits<T>::epsilon()));
 
    const Vec3T<T> delta = a_point - m_tip;
    const Vec2T<T> q(std::sqrt(delta[0] * delta[0] + delta[2] * delta[2]), -delta[1]);
 
    const T d1 = length(q - m_c * Math::max(dot(q, m_c), T(0.0)));
    const T d2 = d1 * ((q[0] * m_c[1] - q[1] * m_c[0] < T(0.0)) ? T(-1.0) : T(1.0));
 
    return d2;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    Vec3T<T> hi = Vec3T<T>::max();
 
    hi[1] = m_tip[1];
 
    return BoundingVolumes::AABBT<T>(Vec3T<T>::min(), hi);
  }
 
protected:
  Vec3T<T> m_tip = Vec3T<T>::zeros();
 
  Vec2T<T> m_c = Vec2T<T>(std::sin(pi<T> / T(8)), std::cos(pi<T> / T(8)));
};
 
template <class T>
class ConeSDF
{
  static_assert(std::is_floating_point_v<T>, "ConeSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  ConeSDF() = default;
 
  ConeSDF(const Vec3T<T>& a_tip, const T& a_height, const T& a_angle) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_tip[0]) && std::isfinite(a_tip[1]) && std::isfinite(a_tip[2]),
                       "ConeSDF: the tip must be finite (%g, %g, %g)",
                       double(a_tip[0]),
                       double(a_tip[1]),
                       double(a_tip[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_height) && a_height > T(0),
                       "ConeSDF: the height must be finite and positive (%g)",
                       double(a_height));
    EBGEOMETRY_REQUIRE(
      a_angle > T(0) && a_angle < T(180), "ConeSDF: the angle must be in (0, 180) degrees (%g)", double(a_angle));
 
    m_tip    = a_tip;
    m_height = a_height;
    m_c[0]   = std::sin(T(0.5) * a_angle * pi<T> / T(180));
    m_c[1]   = std::cos(T(0.5) * a_angle * pi<T> / T(180));
  }
 
  ConeSDF(const ConeSDF&) = default;
 
  ConeSDF(ConeSDF&&) = default;
 
  ConeSDF&
  operator=(const ConeSDF&) = default;
 
  ConeSDF&
  operator=(ConeSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
    EBGEOMETRY_EXPECT(std::abs(length(m_c) - T(1)) < std::sqrt(Math::Limits<T>::epsilon()));
    EBGEOMETRY_EXPECT(m_c[1] > T(0));
 
    const Vec3T<T> delta = a_point - m_tip;
    const T        dr    = std::sqrt(delta[0] * delta[0] + delta[2] * delta[2]);
    const T        dz    = delta[1];
 
    constexpr T zero = T(0.0);
    constexpr T one  = T(1.0);
 
    const Vec2T<T> q = m_height * Vec2T<T>(m_c[0] / m_c[1], -1.0);
    const Vec2T<T> w = Vec2T<T>(dr, dz);
    const Vec2T<T> a = w - Math::clamp(dot(w, q) / dot(q, q), zero, one) * q;
    const Vec2T<T> b = w - Vec2T<T>(q[0] * Math::clamp(w[0] / q[0], zero, one), q[1]);
 
    auto sign = [](const T& x) -> int { return (x > zero) - (x < zero); };
 
    const T k = sign(q[1]);
    const T d = Math::min(dot(a, a), dot(b, b));
    const T s = Math::max(k * (w[0] * q[1] - w[1] * q[0]), k * (w[1] - q[1]));
 
    return std::sqrt(d) * sign(s);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    const T baseRadius = m_height * m_c[0] / m_c[1];
 
    return BoundingVolumes::AABBT<T>(m_tip - Vec3T<T>(baseRadius, m_height, baseRadius),
                                     m_tip + Vec3T<T>(baseRadius, T(0), baseRadius));
  }
 
protected:
  Vec3T<T> m_tip = Vec3T<T>::zeros();
 
  Vec2T<T> m_c = Vec2T<T>(std::sin(pi<T> / T(8)), std::cos(pi<T> / T(8)));
 
  T m_height = T(1);
};
 
template <class T>
class RoundedBoxSDF
{
  static_assert(std::is_floating_point_v<T>, "RoundedBoxSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  RoundedBoxSDF() = default;
 
  RoundedBoxSDF(const Vec3T<T>& a_center, const Vec3T<T>& a_size, const T a_curvature) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center[0]) && std::isfinite(a_center[1]) && std::isfinite(a_center[2]),
                       "RoundedBoxSDF: the center must be finite (%g, %g, %g)",
                       double(a_center[0]),
                       double(a_center[1]),
                       double(a_center[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_curvature) && a_curvature > T(0),
                       "RoundedBoxSDF: the curvature must be finite and positive (%g)",
                       double(a_curvature));
    EBGEOMETRY_REQUIRE(std::isfinite(a_size[0]) && std::isfinite(a_size[1]) && std::isfinite(a_size[2]) &&
                         a_size[0] > T(2) * a_curvature && a_size[1] > T(2) * a_curvature &&
                         a_size[2] > T(2) * a_curvature,
                       "RoundedBoxSDF: the size (%g, %g, %g) must be finite and exceed twice the curvature (%g)",
                       double(a_size[0]),
                       double(a_size[1]),
                       double(a_size[2]),
                       double(a_curvature));
 
    m_center     = a_center;
    m_dimensions = T(0.5) * a_size - a_curvature * Vec3T<T>::ones();
    m_curvature  = a_curvature;
  }
 
  RoundedBoxSDF(const RoundedBoxSDF&) = default;
 
  RoundedBoxSDF(RoundedBoxSDF&&) = default;
 
  RoundedBoxSDF&
  operator=(const RoundedBoxSDF&) = default;
 
  RoundedBoxSDF&
  operator=(RoundedBoxSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    // q is the point's per-axis distance beyond the inner box (negative inside it). Outside, the
    // distance to the inner box is the length of q's positive part; inside, it is the distance to
    // the nearest face, max(q) <= 0. Either way the rounding moves the surface out by m_curvature.
    const Vec3T<T> p = a_point - m_center;
    const Vec3T<T> q(
      std::abs(p[0]) - m_dimensions[0], std::abs(p[1]) - m_dimensions[1], std::abs(p[2]) - m_dimensions[2]);
 
    const T outside = length(max(q, Vec3T<T>::zeros()));
    const T inside  = Math::min(Math::max(q[0], Math::max(q[1], q[2])), T(0));
 
    return outside + inside - m_curvature;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    const Vec3T<T> half = m_dimensions + m_curvature * Vec3T<T>::ones();
 
    return BoundingVolumes::AABBT<T>(m_center - half, m_center + half);
  }
 
protected:
  T m_curvature = T(0.1);
 
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  Vec3T<T> m_dimensions = Vec3T<T>(T(0.4), T(0.4), T(0.4));
};
 
template <class T>
class PerlinSDF
{
  static_assert(std::is_floating_point_v<T>, "PerlinSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::NotADistance;
 
  PerlinSDF() noexcept : PerlinSDF(T(1), Vec3T<T>::ones(), T(0.5), 1U)
  {}
 
  PerlinSDF(const T            a_noiseAmplitude,
            const Vec3T<T>     a_noiseFrequency,
            const T            a_noisePersistence,
            const unsigned int a_noiseOctaves) noexcept
  {
    EBGEOMETRY_REQUIRE(
      std::isfinite(a_noiseAmplitude), "PerlinSDF: the noise amplitude must be finite (%g)", double(a_noiseAmplitude));
    EBGEOMETRY_REQUIRE(std::isfinite(a_noiseFrequency[0]) && std::isfinite(a_noiseFrequency[1]) &&
                         std::isfinite(a_noiseFrequency[2]),
                       "PerlinSDF: the noise frequency must be finite (%g, %g, %g)",
                       double(a_noiseFrequency[0]),
                       double(a_noiseFrequency[1]),
                       double(a_noiseFrequency[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_noisePersistence),
                       "PerlinSDF: the noise persistence must be finite (%g)",
                       double(a_noisePersistence));
 
    m_noiseAmplitude   = a_noiseAmplitude;
    m_noiseFrequency   = a_noiseFrequency;
    m_noisePersistence = Math::min(T(1), a_noisePersistence);
    m_noiseOctaves     = Math::max(1U, a_noiseOctaves);
 
    // Each octave divides the frequency by the persistence, so zero would give an infinite frequency
    // (and NaN noise) for an octave that contributes nothing anyway.
    if (!(m_noisePersistence > T(0))) {
      m_noisePersistence = T(1);
      m_noiseOctaves     = 1U;
    }
 
    std::mt19937 g(0);
    this->shuffle(g);
  }
 
  PerlinSDF(const PerlinSDF&) = default;
 
  PerlinSDF(PerlinSDF&&) = default;
 
  PerlinSDF&
  operator=(const PerlinSDF&) = default;
 
  PerlinSDF&
  operator=(PerlinSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    T ret = 0.0;
 
    Vec3T<T> curFreq = m_noiseFrequency;
    T        curAmp  = 1.0;
 
    for (unsigned int curOctave = 0; curOctave < m_noiseOctaves; curOctave++) {
      ret += T(0.5) * curAmp * (T(1) + this->noise(a_point * curFreq));
 
      curFreq = curFreq / m_noisePersistence;
      curAmp  = curAmp * m_noisePersistence;
    }
 
    return ret * m_noiseAmplitude;
  };
 
  template <class URNG>
  void
  shuffle(URNG& g) noexcept
  {
 
    for (unsigned int i = 0; i < 256; i++) {
      m_permutationTable[i] = static_cast<int>(i);
    }
 
    std::shuffle(m_permutationTable.begin(), m_permutationTable.begin() + 256, g);
 
    for (int i = 0; i < 256; i++) {
      m_permutationTable[i + 256] = m_permutationTable[i];
    }
  }
 
  Array<int, 512>&
  getPermutationTable() noexcept
  {
    return m_permutationTable;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    return BoundingVolumes::AABBT<T>(Vec3T<T>::min(), Vec3T<T>::max());
  }
 
protected:
  Vec3T<T> m_noiseFrequency = Vec3T<T>(T(1), T(1), T(1));
 
  T m_noiseAmplitude = T(1);
 
  T m_noisePersistence = T(0.5);
 
  Array<int, 512> m_permutationTable = {};
 
  unsigned int m_noiseOctaves = 1U;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  lerp(const T t, const T a, const T b) const noexcept
  {
    return a + t * (b - a);
  };
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  fade(const T t) const noexcept
  {
    return t * t * t * (t * (t * 6 - 15) + 10);
  };
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  grad(const int hash, const T x, const T y, const T z) const noexcept
  {
    const int h = hash & 15;
    const T   u = h < 8 ? x : y;
    const T   v = h < 4 ? y : h == 12 || h == 14 ? x : z;
    return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  noise(const Vec3T<T>& a_point) const noexcept
  {
    // Lower cube corner
    const int X = static_cast<int>(std::floor(a_point[0])) & 255;
    const int Y = static_cast<int>(std::floor(a_point[1])) & 255;
    const int Z = static_cast<int>(std::floor(a_point[2])) & 255;
 
    // Relative distance wrt lower cube corner
    const T x = a_point[0] - std::floor(a_point[0]);
    const T y = a_point[1] - std::floor(a_point[1]);
    const T z = a_point[2] - std::floor(a_point[2]);
 
    // Fade curves
    const T u = fade(x);
    const T v = fade(y);
    const T w = fade(z);
 
    // Hash coordinates of 8 cube corners
    const int A  = m_permutationTable[X] + Y;
    const int AA = m_permutationTable[A] + Z;
    const int AB = m_permutationTable[A + 1] + Z;
    const int B  = m_permutationTable[X + 1] + Y;
    const int BA = m_permutationTable[B] + Z;
    const int BB = m_permutationTable[B + 1] + Z;
 
    // Add blended results from 8 corners of cube
    return lerp(
      w,
      lerp(v,
           lerp(u, grad(m_permutationTable[AA], x, y, z), grad(m_permutationTable[BA], x - 1, y, z)),
           lerp(u, grad(m_permutationTable[AB], x, y - 1, z), grad(m_permutationTable[BB], x - 1, y - 1, z))),
      lerp(v,
           lerp(u, grad(m_permutationTable[AA + 1], x, y, z - 1), grad(m_permutationTable[BA + 1], x - 1, y, z - 1)),
           lerp(u,
                grad(m_permutationTable[AB + 1], x, y - 1, z - 1),
                grad(m_permutationTable[BB + 1], x - 1, y - 1, z - 1))));
  };
};
 
template <class T>
class RoundedCylinderSDF
{
  static_assert(std::is_floating_point_v<T>, "RoundedCylinderSDF<T>: T must be a floating-point type");
 
public:
  static constexpr DistanceQuality distanceQuality = DistanceQuality::Exact;
 
  RoundedCylinderSDF() = default;
 
  RoundedCylinderSDF(const Vec3T<T>& a_center, const T a_radius, const T a_curvature, const T a_height) noexcept
  {
    EBGEOMETRY_REQUIRE(std::isfinite(a_center[0]) && std::isfinite(a_center[1]) && std::isfinite(a_center[2]),
                       "RoundedCylinderSDF: the center must be finite (%g, %g, %g)",
                       double(a_center[0]),
                       double(a_center[1]),
                       double(a_center[2]));
    EBGEOMETRY_REQUIRE(std::isfinite(a_radius) && a_radius > T(0),
                       "RoundedCylinderSDF: the radius must be finite and positive (%g)",
                       double(a_radius));
    EBGEOMETRY_REQUIRE(std::isfinite(a_curvature) && a_curvature > T(0),
                       "RoundedCylinderSDF: the curvature must be finite and positive (%g)",
                       double(a_curvature));
    EBGEOMETRY_REQUIRE(std::isfinite(a_height) && a_height > T(0),
                       "RoundedCylinderSDF: the height must be finite and positive (%g)",
                       double(a_height));
    EBGEOMETRY_REQUIRE(a_curvature < a_radius,
                       "RoundedCylinderSDF: the curvature (%g) must be less than the radius (%g)",
                       double(a_curvature),
                       double(a_radius));
    EBGEOMETRY_REQUIRE(T(2) * a_curvature < a_height,
                       "RoundedCylinderSDF: twice the curvature (%g) must be less than the height (%g)",
                       double(a_curvature),
                       double(a_height));
 
    m_center      = a_center;
    m_majorRadius = a_radius - a_curvature;
    m_minorRadius = a_curvature;
    m_height      = T(0.5) * a_height - a_curvature;
  }
 
  RoundedCylinderSDF(const RoundedCylinderSDF&) = default;
 
  RoundedCylinderSDF(RoundedCylinderSDF&&) = default;
 
  RoundedCylinderSDF&
  operator=(const RoundedCylinderSDF&) = default;
 
  RoundedCylinderSDF&
  operator=(RoundedCylinderSDF&&) = default;
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  T
  signedDistance(const Vec3T<T>& a_point) const noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    const Vec3T<T> p  = a_point - m_center;
    const T        xz = std::sqrt(p[2] * p[2] + p[0] * p[0]);
    const auto     d1 = Vec2T<T>(xz - m_majorRadius, std::abs(p[1]) - m_height);
    const auto     d2 = Vec2T<T>(Math::max(d1[0], T(0)), Math::max(d1[1], T(0)));
 
    return Math::min(Math::max(d1[0], d1[1]), T(0)) + std::sqrt(d2[0] * d2[0] + d2[1] * d2[1]) - m_minorRadius;
  }
 
  [[nodiscard]] EBGEOMETRY_HOST_DEVICE
  BoundingVolumes::AABBT<T>
  computeBoundingVolume() const noexcept
  {
    const T        radius = m_majorRadius + m_minorRadius;
    const Vec3T<T> half(radius, m_height + m_minorRadius, radius);
 
    return BoundingVolumes::AABBT<T>(m_center - half, m_center + half);
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_majorRadius = T(0.9);
 
  T m_minorRadius = T(0.1);
 
  T m_height = T(0.4);
};
 
static_assert(std::is_trivially_copyable_v<PlaneSDF<float>>, "PlaneSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<PlaneSDF<double>>, "PlaneSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<SphereSDF<float>>, "SphereSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<SphereSDF<double>>, "SphereSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<BoxSDF<float>>, "BoxSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<BoxSDF<double>>, "BoxSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<TorusSDF<float>>, "TorusSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<TorusSDF<double>>, "TorusSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<CylinderSDF<float>>, "CylinderSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<CylinderSDF<double>>, "CylinderSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<InfiniteCylinderSDF<float>>,
              "InfiniteCylinderSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<InfiniteCylinderSDF<double>>,
              "InfiniteCylinderSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<CapsuleSDF<float>>, "CapsuleSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<CapsuleSDF<double>>, "CapsuleSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<InfiniteConeSDF<float>>,
              "InfiniteConeSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<InfiniteConeSDF<double>>,
              "InfiniteConeSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<ConeSDF<float>>, "ConeSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<ConeSDF<double>>, "ConeSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<RoundedBoxSDF<float>>, "RoundedBoxSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<RoundedBoxSDF<double>>, "RoundedBoxSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<PerlinSDF<float>>, "PerlinSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<PerlinSDF<double>>, "PerlinSDF<double> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<RoundedCylinderSDF<float>>,
              "RoundedCylinderSDF<float> must be trivially copyable");
static_assert(std::is_trivially_copyable_v<RoundedCylinderSDF<double>>,
              "RoundedCylinderSDF<double> must be trivially copyable");
 
} // namespace EBGeometry
 
// Every translation unit that sees these shapes sees their ray forms too, so firstHit() on one of them
// is the same function everywhere.
#include "EBGeometry_AnalyticRayForms.hpp"
 
#endif