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 <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>
#include <random>
#include <type_traits>
 
// Our includes
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_Constants.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_SignedDistanceFunction.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
template <class T>
class PlaneSDF : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "PlaneSDF<T>: T must be a floating-point type");
 
public:
  PlaneSDF() = default;
 
  PlaneSDF(const Vec3T<T>& a_point, const Vec3T<T>& a_normal) 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::isfinite(a_normal[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_normal[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_normal[2]));
    EBGEOMETRY_EXPECT(a_normal.length() > T(0));
 
    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;
 
  ~PlaneSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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(std::numeric_limits<T>::epsilon()));
 
    return dot((a_point - m_point), 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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "SphereSDF<T>: T must be a floating-point type");
 
public:
  SphereSDF() = default;
 
  SphereSDF(const Vec3T<T>& a_center, const T& a_radius) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_center[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_radius));
    EBGEOMETRY_EXPECT(a_radius > T(0));
 
    m_center = a_center;
    m_radius = a_radius;
  }
 
  SphereSDF(const SphereSDF&) = default;
 
  SphereSDF(SphereSDF&&) = default;
 
  SphereSDF&
  operator=(const SphereSDF&) = default;
 
  SphereSDF&
  operator=(SphereSDF&&) = default;
 
  ~SphereSDF() override = default;
 
  [[nodiscard]] const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  [[nodiscard]] const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  T&
  getRadius() noexcept
  {
    return m_radius;
  }
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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;
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_radius = T(1);
};
 
template <class T>
class BoxSDF : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "BoxSDF<T>: T must be a floating-point type");
 
public:
  BoxSDF() = default;
 
  BoxSDF(const Vec3T<T>& a_loCorner, const Vec3T<T>& a_hiCorner) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_loCorner[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_loCorner[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_loCorner[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_hiCorner[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_hiCorner[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_hiCorner[2]));
    EBGEOMETRY_EXPECT(a_loCorner[0] < a_hiCorner[0]);
    EBGEOMETRY_EXPECT(a_loCorner[1] < a_hiCorner[1]);
    EBGEOMETRY_EXPECT(a_loCorner[2] < 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;
 
  ~BoxSDF() override = default;
 
  [[nodiscard]] const Vec3T<T>&
  getLowCorner() const noexcept
  {
    return m_loCorner;
  }
 
  Vec3T<T>&
  getLowCorner() noexcept
  {
    return m_loCorner;
  }
 
  [[nodiscard]] const Vec3T<T>&
  getHighCorner() const noexcept
  {
    return m_hiCorner;
  }
 
  Vec3T<T>&
  getHighCorner() noexcept
  {
    return m_hiCorner;
  }
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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(std::max(m_loCorner[0] - a_point[0], a_point[0] - m_hiCorner[0]),
                         std::max(m_loCorner[1] - a_point[1], a_point[1] - m_hiCorner[1]),
                         std::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
    // std::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 = std::min(T(0.0), delta[delta.maxDir(false)]) + max(Vec3T<T>::zeros(), delta).length();
 
    return d;
  }
 
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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "TorusSDF<T>: T must be a floating-point type");
 
public:
  TorusSDF() = default;
 
  TorusSDF(const Vec3T<T>& a_center, const T& a_majorRadius, const T& a_minorRadius) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_center[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_majorRadius));
    EBGEOMETRY_EXPECT(std::isfinite(a_minorRadius));
    EBGEOMETRY_EXPECT(a_majorRadius > T(0));
    EBGEOMETRY_EXPECT(a_minorRadius > T(0));
    EBGEOMETRY_EXPECT(a_minorRadius < 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;
 
  ~TorusSDF() override = default;
 
  [[nodiscard]] const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  [[nodiscard]] const T&
  getMajorRadius() const noexcept
  {
    return m_majorRadius;
  }
 
  T&
  getMajorRadius() noexcept
  {
    return m_majorRadius;
  }
 
  [[nodiscard]] const T&
  getMinorRadius() const noexcept
  {
    return m_minorRadius;
  }
 
  T&
  getMinorRadius() noexcept
  {
    return m_minorRadius;
  }
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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[1] * p[1]) - m_majorRadius;
    const T        d   = std::sqrt(rho * rho + p[2] * p[2]) - m_minorRadius;
 
    return d;
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_majorRadius = T(1);
 
  T m_minorRadius = T(0.5);
};
 
template <class T>
class CylinderSDF : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "CylinderSDF<T>: T must be a floating-point type");
 
public:
  CylinderSDF() = default;
 
  CylinderSDF(const Vec3T<T>& a_center1, const Vec3T<T>& a_center2, const T& a_radius) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_center1[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center1[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center1[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center2[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center2[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center2[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_radius));
    EBGEOMETRY_EXPECT(a_radius > T(0));
    EBGEOMETRY_EXPECT((a_center2 - a_center1).length() > T(0));
 
    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;
 
  ~CylinderSDF() override = default;
 
  [[nodiscard]] const Vec3T<T>&
  getCenter1() const noexcept
  {
    return m_center1;
  }
 
  [[nodiscard]] const Vec3T<T>&
  getCenter2() const noexcept
  {
    return m_center2;
  }
 
  [[nodiscard]] const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    T d = std::numeric_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(std::numeric_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;
  }
 
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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "InfiniteCylinderSDF<T>: T must be a floating-point type");
 
public:
  InfiniteCylinderSDF() = default;
 
  InfiniteCylinderSDF(const Vec3T<T>& a_center, const T& a_radius, const size_t a_axis) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_center[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_center[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_radius));
    EBGEOMETRY_EXPECT(a_radius > T(0));
    EBGEOMETRY_EXPECT(a_axis < 3U);
 
    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;
 
  ~InfiniteCylinderSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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;
  }
 
protected:
  Vec3T<T> m_center = Vec3T<T>::zeros();
 
  T m_radius = T(1);
 
  size_t m_axis = 1U;
};
 
template <class T>
class CapsuleSDF : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "CapsuleSDF<T>: T must be a floating-point type");
 
public:
  CapsuleSDF() = default;
 
  CapsuleSDF(const Vec3T<T>& a_tip1, const Vec3T<T>& a_tip2, const T& a_radius) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_tip1[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip1[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip1[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip2[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip2[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip2[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_radius));
    EBGEOMETRY_EXPECT(a_radius > T(0));
    EBGEOMETRY_EXPECT((a_tip2 - a_tip1).length() > T(0));
 
    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;
 
  ~CapsuleSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
    EBGEOMETRY_EXPECT((m_center2 - m_center1).length() > T(0));
 
    const Vec3T<T> v1 = a_point - m_center1;
    const Vec3T<T> v2 = m_center2 - m_center1;
 
    const T h = std::clamp(dot(v1, v2) / dot(v2, v2), T(0.0), T(1.0));
    const T d = length(v1 - h * v2) - m_radius;
 
    return d;
  }
 
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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "InfiniteConeSDF<T>: T must be a floating-point type");
 
public:
  InfiniteConeSDF() = default;
 
  InfiniteConeSDF(const Vec3T<T>& a_tip, const T& a_angle) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_angle));
    EBGEOMETRY_EXPECT(a_angle > T(0));
    EBGEOMETRY_EXPECT(a_angle < T(180));
 
    m_tip = a_tip;
    m_c.x = std::sin(T(0.5) * a_angle * pi<T> / T(180));
    m_c.y = 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;
 
  ~InfiniteConeSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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(std::numeric_limits<T>::epsilon()));
 
    const Vec3T<T> delta = a_point - m_tip;
    const Vec2T<T> q(std::sqrt(delta[0] * delta[0] + delta[1] * delta[1]), -delta[2]);
 
    const T d1 = length(q - m_c * std::max(dot(q, m_c), T(0.0)));
    const T d2 = d1 * ((q.x * m_c.y - q.y * m_c.x < T(0.0)) ? T(-1.0) : T(1.0));
 
    return d2;
  }
 
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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "ConeSDF<T>: T must be a floating-point type");
 
public:
  ConeSDF() = default;
 
  ConeSDF(const Vec3T<T>& a_tip, const T& a_height, const T& a_angle) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_tip[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_height));
    EBGEOMETRY_EXPECT(std::isfinite(a_angle));
    EBGEOMETRY_EXPECT(a_height > T(0));
    EBGEOMETRY_EXPECT(a_angle > T(0));
    EBGEOMETRY_EXPECT(a_angle < T(180));
 
    m_tip    = a_tip;
    m_height = a_height;
    m_c.x    = std::sin(T(0.5) * a_angle * pi<T> / T(180));
    m_c.y    = 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;
 
  ~ConeSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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(std::numeric_limits<T>::epsilon()));
    EBGEOMETRY_EXPECT(m_c.y > T(0));
 
    const Vec3T<T> delta = a_point - m_tip;
    const T        dr    = std::sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
    const T        dz    = delta[2];
 
    constexpr T zero = T(0.0);
    constexpr T one  = T(1.0);
 
    const Vec2T<T> q = m_height * Vec2T<T>(m_c.x / m_c.y, -1.0);
    const Vec2T<T> w = Vec2T<T>(dr, dz);
    const Vec2T<T> a = w - std::clamp(dot(w, q) / dot(q, q), zero, one) * q;
    const Vec2T<T> b = w - Vec2T<T>(q.x * std::clamp(w.x / q.x, zero, one), q.y);
 
    auto sign = [](const T& x) -> int { return (x > zero) - (x < zero); };
 
    const T k = sign(q.y);
    const T d = std::min(dot(a, a), dot(b, b));
    const T s = std::max(k * (w.x * q.y - w.y * q.x), k * (w.y - q.y));
 
    return std::sqrt(d) * sign(s);
  }
 
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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "RoundedBoxSDF<T>: T must be a floating-point type");
 
public:
  RoundedBoxSDF() = default;
 
  RoundedBoxSDF(const Vec3T<T>& a_dimensions, const T a_curvature) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_dimensions[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_dimensions[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_dimensions[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_curvature));
    EBGEOMETRY_EXPECT(a_dimensions[0] > T(0));
    EBGEOMETRY_EXPECT(a_dimensions[1] > T(0));
    EBGEOMETRY_EXPECT(a_dimensions[2] > T(0));
    EBGEOMETRY_EXPECT(a_curvature > T(0));
 
    m_dimensions = T(0.5) * a_dimensions;
    m_sphere     = std::make_shared<SphereSDF<T>>(Vec3T<T>::zeros(), a_curvature);
  }
 
  RoundedBoxSDF(const RoundedBoxSDF&) = default;
 
  RoundedBoxSDF(RoundedBoxSDF&&) = default;
 
  RoundedBoxSDF&
  operator=(const RoundedBoxSDF&) = default;
 
  RoundedBoxSDF&
  operator=(RoundedBoxSDF&&) = default;
 
  ~RoundedBoxSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    return m_sphere->signedDistance(a_point - clamp(a_point, -m_dimensions, m_dimensions));
  }
 
protected:
  std::shared_ptr<SphereSDF<T>> m_sphere = std::make_shared<SphereSDF<T>>(Vec3T<T>::zeros(), T(0.1));
 
  Vec3T<T> m_dimensions = Vec3T<T>(T(0.5), T(0.5), T(0.5));
};
 
template <class T>
class PerlinSDF : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "PerlinSDF<T>: T must be a floating-point type");
 
public:
  PerlinSDF() = default;
 
  PerlinSDF(const T            a_noiseAmplitude,
            const Vec3T<T>     a_noiseFrequency,
            const T            a_noisePersistence,
            const unsigned int a_noiseOctaves) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_noiseAmplitude));
    EBGEOMETRY_EXPECT(std::isfinite(a_noiseFrequency[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_noiseFrequency[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_noiseFrequency[2]));
    EBGEOMETRY_EXPECT(std::isfinite(a_noisePersistence));
 
    m_noiseAmplitude   = a_noiseAmplitude;
    m_noiseFrequency   = a_noiseFrequency;
    m_noisePersistence = std::min(T(1), a_noisePersistence);
    m_noiseOctaves     = std::max(1U, a_noiseOctaves);
 
    for (int i = 0; i < 256; i++) {
      m_permutationTable[i]       = s_perlinPermutationTable[i];
      m_permutationTable[i + 256] = s_perlinPermutationTable[i];
    }
 
    std::mt19937 g(0);
    this->shuffle(g);
  }
 
  PerlinSDF(const PerlinSDF&) = default;
 
  PerlinSDF(PerlinSDF&&) = default;
 
  PerlinSDF&
  operator=(const PerlinSDF&) = default;
 
  PerlinSDF&
  operator=(PerlinSDF&&) = default;
 
  ~PerlinSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    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];
    }
  }
 
  std::array<int, 512>&
  getPermutationTable() noexcept
  {
    return m_permutationTable;
  }
 
protected:
  static constexpr std::array<int, 256> s_perlinPermutationTable = {
    151, 160, 137, 91,  90,  15,  131, 13,  201, 95,  96,  53,  194, 233, 7,   225, 140, 36,  103, 30,  69,  142,
    8,   99,  37,  240, 21,  10,  23,  190, 6,   148, 247, 120, 234, 75,  0,   26,  197, 62,  94,  252, 219, 203,
    117, 35,  11,  32,  57,  177, 33,  88,  237, 149, 56,  87,  174, 20,  125, 136, 171, 168, 68,  175, 74,  165,
    71,  134, 139, 48,  27,  166, 77,  146, 158, 231, 83,  111, 229, 122, 60,  211, 133, 230, 220, 105, 92,  41,
    55,  46,  245, 40,  244, 102, 143, 54,  65,  25,  63,  161, 1,   216, 80,  73,  209, 76,  132, 187, 208, 89,
    18,  169, 200, 196, 135, 130, 116, 188, 159, 86,  164, 100, 109, 198, 173, 186, 3,   64,  52,  217, 226, 250,
    124, 123, 5,   202, 38,  147, 118, 126, 255, 82,  85,  212, 207, 206, 59,  227, 47,  16,  58,  17,  182, 189,
    28,  42,  223, 183, 170, 213, 119, 248, 152, 2,   44,  154, 163, 70,  221, 153, 101, 155, 167, 43,  172, 9,
    129, 22,  39,  253, 19,  98,  108, 110, 79,  113, 224, 232, 178, 185, 112, 104, 218, 246, 97,  228, 251, 34,
    242, 193, 238, 210, 144, 12,  191, 179, 162, 241, 81,  51,  145, 235, 249, 14,  239, 107, 49,  192, 214, 31,
    181, 199, 106, 157, 184, 84,  204, 176, 115, 121, 50,  45,  127, 4,   150, 254, 138, 236, 205, 93,  222, 114,
    67,  29,  24,  72,  243, 141, 128, 195, 78,  66,  215, 61,  156, 180};
 
  Vec3T<T> m_noiseFrequency = Vec3T<T>(T(1), T(1), T(1));
 
  T m_noiseAmplitude = T(1);
 
  T m_noisePersistence = T(0.5);
 
  std::array<int, 512> m_permutationTable = {};
 
  unsigned int m_noiseOctaves = 1U;
 
  [[nodiscard]] virtual T
  lerp(const T t, const T a, const T b) const noexcept
  {
    return a + t * (b - a);
  };
 
  [[nodiscard]] virtual T
  fade(const T t) const noexcept
  {
    return t * t * t * (t * (t * 6 - 15) + 10);
  };
 
  [[nodiscard]] 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]] 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 : public SignedDistanceFunction<T>
{
  static_assert(std::is_floating_point_v<T>, "RoundedCylinderSDF<T>: T must be a floating-point type");
 
public:
  RoundedCylinderSDF() = default;
 
  RoundedCylinderSDF(const T a_radius, const T a_curvature, const T a_height) noexcept
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_radius));
    EBGEOMETRY_EXPECT(std::isfinite(a_curvature));
    EBGEOMETRY_EXPECT(std::isfinite(a_height));
    EBGEOMETRY_EXPECT(a_radius > T(0));
    EBGEOMETRY_EXPECT(a_curvature > T(0));
    EBGEOMETRY_EXPECT(a_height > T(0));
    EBGEOMETRY_EXPECT(a_curvature < a_radius);
    EBGEOMETRY_EXPECT(T(2) * a_curvature < a_height);
 
    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;
 
  ~RoundedCylinderSDF() override = default;
 
  [[nodiscard]] T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    EBGEOMETRY_EXPECT(std::isfinite(a_point[0]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[1]));
    EBGEOMETRY_EXPECT(std::isfinite(a_point[2]));
 
    const T    xz = std::sqrt(a_point[2] * a_point[2] + a_point[0] * a_point[0]);
    const auto d1 = Vec2T<T>(xz - m_majorRadius, std::abs(a_point[1]) - m_height);
    const auto d2 = Vec2T<T>(std::max(d1.x, T(0)), std::max(d1.y, T(0)));
 
    return std::min(std::max(d1.x, d1.y), T(0)) + std::sqrt(d2.x * d2.x + d2.y * d2.y) - m_minorRadius;
  }
 
protected:
  T m_majorRadius = T(0.9);
 
  T m_minorRadius = T(0.1);
 
  T m_height = T(0.4);
};
 
} // namespace EBGeometry
 
#endif